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ABSTRACT 


This  report  represents  a  study  of  the 
equations  governing  the  propagation  of  weakly  nonlinear 
waves  in  regions  where  currents  exist  and  where  the  depth 
and  current  are  allowed  to  vary.  After  deriving  a  velocity 
potential  applicable  to  propagating  surface  waves  according 
to  the  Stoxes  expansion,  the  Lagrangian  governing  wave 
motion  in  the  propagation  space  is  derived.  A  consistent 
perturbation  scheme  then  leads  to  the  equations  governing 
the  wave  and  current  motion  at  each  order  in  the  Stokes 
series . 

After  neglecting  time  dependence  of  the  wave 
amplitude,  parabolic  equations  governing  the  combined 
refraction  and  diffraction  of  Stokes  waves  of  small 
amplitude  are  developed  and  used  to  calculate  the  wavefields 
for  several  representative  cases  illustrating  the  importance 
of  nonlinear  effects.  The  computational  models  are  verified 
by  comparison  to  laboratory  data  of  wave  amplitude  for  a 
wavefield  focussed  by  a  submerged  shoal,  and  it  is  found 


xi 


that  the  nonlinear  model  accounts  for  the  nrsa^or 
discrepancies  found  between  linear  model  results  and  tne 
laboratory  data. 


tl 

ii 


i 

Mil 


CHAPTER  1.  INTRODUCTION 


One  of  the  central  problems  in  the  study  of  the 
coastal  surface  wave  environment  is  the  prediction  of  the 
transformation  of  waves  as  they  propagate  from  the  deep 
ocean  towards  tne  shore.  During  this  period  of  propagation, 
the  waves  are  subject  to  many  of  the  physical  effects 
characteristic  of  general  wave  problems  in  many  branches  of 
physics.  First,  the  domain  of  propagation  may  in  general  be 
in  motion,  due  to  the  effect  of  tidal  or  other  currents. 
For  the  case  of  currents  with  scales  of  variation  much 
larger  than  the  scale  of  the  local  wave,  which  may  be 
characterised  by  the  wavelength,  this  effect  may  be 
quantified  locally  in  terms  of  a  simple  Galilean 
transformation  for  waves  of  small  amplitude,  in  which  the 
frequency  of  the  wave  oscillation  is  defined  with  respect  to 
the  moving  medium.  Secondly,  the  propagation  medium  is,  in 
general,  inhomogeneous.  The  innomogenei ty  may  arise  locally 
due  to  a  smaller  scale  of  variation  of  the  ambient  current, 
and  will  be  apparent  over  larger  distances  for  any  scale  of 
variation.  Additionally,  after  the  waves  propagate  into 


1 


2 


water  snallow  enough  so  that  tne  phase  velocity  of  the  wave 
becomes  a  function  of  the  local  depth,  the  pattern  of  wave 
propagation  will  be  affected  by  the  spatial  variations  of 
depth.  Each  of  these  inhomogeneities  can  cause  refractive 
and  diffractive  effects,  leading  to  the  presence  of  caustics 
and  related  characteristics  of  the  wave  field. 

In  addition  to  effects  related  to  movement  and 
spatial  inhomogeneity  of  the  domain,  the  wave  motion  itself 
is  nonlinear,  and  thus  the  components  of  the  wave  field  may 
interact  in  both  subtle  and  spectacular  ways.  The  simplest 
interaction  involves  the  mutual  interaction  of  the 
components  of  a  single  wave,  discussed  by  Phillips  (1960), 
leading  to  the  effect  of  amplitude  dispersion,  in  which  tne 
phase  speed  of  the  wave  is  increased  as  a  function  of  its 
local  amplitude.  In  addition,  separate  wave  trains,  either 
propagating  colinearly  or  intersecting  obliquely,  may 
exchange  energy  in  a  complex  manner  leading  to  gradual 
shifts  of  the  overall  characteristics  of  the  wave  field. 
The  individual  waves  of  the  wave  field  may  respond  to  the 
interactions  in  an  unstable  manner,  leading  to  amplitude 
modulations  which  act  over  a  smail  spatial  scale  and  which 
also  effect  the  phase  velocity  of  the  individual  components. 

Wave  and  current  motions  are  subject  both  to  direct 
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viscous  effects  and  to  small  scale  turbulence  which  may  be 
modelled  in  analogy  to  viscosity.  In  general,  the  presence 
of  viscous  effects  causes  the  domain  of  propagation  to  be 
dissipative,  so  that  waves  propagating  over  long  distances 
may  lose  a  significant  proportion  of  their  energy  content, 
principally  through  interaction  with  the  bottom,  viscosity, 
as  well  as  the  rotation  of  the  earth,  also  leads  to 
rotational  current  distributions.  Finally,  the  wavefield  is 
subject  to  the  action  of  applied  surface  shear  stress  due  to 
wind,  and  to  spatial  fluctuations  in  the  applied  atmospheric 
pressure.  The  shear  stress  due  to  wind  is  important  in  the 
scale  of  the  relatively  snort  surface  waves,  and  leads  to 
surface  drift  currents  as  well  as  growth  of  the  waves  as 
energy  is  transferred  to  the  water  column. 

In  order  for  an  analytic  model  of  wave  propagation 
to  be  successfully  applicable  over  all  propagation  scales, 
it  must  oegin  to  address  all  of  the  effects  listed  above. 
The  general  problem  of  water  wave  propagation  is  therefore 
extremely  complex,  and  progress  in  its  study  has  necessarily 
been  made  through  the  study  of  subsets  of  the  general 
problem  applicable  to  problems  of  a  more  simplified  nature. 
A  major  simplification  which  still  allows  for  a  general 
treatment  of  a  wide  range  of  the  mentioned  effects  is  the 
neglect  of  viscosity  and  wind  shear  effects,  leading  to  the 
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study  of  wave  states  which  do  not  dissipate  or  increase  in 
energy.  This  simplification  leads  to  the  ability  to  study  a 
subset  of  the  governing  equations  applicable  to 
incompressible  and  inviscid  flows.  The  further  assumption 
of  irrotationality  allows  the  entire  fluid  flow  to  be 
described  by  a  velocity  potential.  It  should  be  noted  that 
the  assumption  of  irrotationality,  especially  with  regard  to 
the  large  scale  current,  is  strictly  invalid  in  most 
physical  cases,  and  it  would  be  desirable  to  obtain  an 
estimate  of  the  error  incurred  by  the  neglect  of  this  effect 
in  any  given  case.  In  this  study,  the  assumption  of 
irrotationality  will  be  included.  The  remaining  subset  of 
the  problem  still  allows  for  the  study  of  propagation  of 
nonlinear  waves  through  an  inhomogeneous,  moving  domain, 
which  is  the  problem  investigated  here.  The  governing 
differential  equations  and  a  corresponding  variational 
principle  for  inviscid  and  irrotational  flow  are  discussed 
in  Chapter  2,  along  with  the  scales  of  variation  to  be 
studied  and  the  small  parameter  expansion  of  the  independent 
variables . 

In  most  instances  to  date,  the  problems  of  nonlinear 
wave-wave  interaction  and  refraction-  diffraction  have  been 
studied  as  separate  entities.  After  linearization  of  the 
free  surface  boundary  conditions,  the  sea  state  may  be 
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characterized  as  the  spectral  sum  of  individual, 
non-interacting  linear  components,  with  the  consequence  that 
the  propagation  characteristics  of  each  component  may  be 
calculated  separately.  As  a  consequence,  much  of  the  work 
on  wave-current  interaction  and  refraction  -diffraction  nas 
centered  on  the  study  of  single  component  wave  systems, 
where  the  time  dependence  may  be  characterised  entirely  in 
terms  of  the  absolute  wave  frequency  within  a  stationary 
reference  frame.  Further,  the  problems  of  refraction  and 
diffraction  have  historically  been  considered  separately  and 
have  been  joined  successfully  only  in  the  past  decade.  The 
first  method  for  studying  surface  wave  refraction  in  an 
arbitrarily  varying  domain  was  proposed  by  Obrien  and 
Mason (1940) ,  and  is  referred  to  as  the  ray  method.  Rays,  or 
ortnogonals  to  the  travelling  wave  crest,  are  traced 
independently  through  the  domain  and  are  assumed  to  be 
everywhere  parallel  to  the  local  direction  of  wave 
propagation.  The  local  wave  amplitude  is  then  calculated 
based  on  the  assumption  that  energy  does  not  cross  the  rays; 
i.e.,  that  the  area  between  adjacent  rays  is  a  tube  of 
constant  energy  flux.  Arthur  (1950)  extended  the  ray  method 
to  include  the  effect  of  wave-current  interaction,  and 
calculated  the  rays  by  means  of  Fermat's  principle  of  least 
time.  Skovgaard,  Jonsson  and  Bertelsen  (1975)  extended  the 
method  to  include  the  effect  of  wave-energy  dissipation  due 


to  bottom  friction. 


Keller(1958)  showed  that  the  refraction 
approximation  arises  as  the  first,  or  geometric  optics, 
approximation  in  a  WKB  expansion  of  the  governing  equations, 
after  neglecting  the  effect  of  wave  amplitude  gradients 
normal  to  the  rays.  The  limitation  of  the  approximation 
became  apparent  in  the  case  where  one  or  more  rays 
intersect,  leading  to  singularities  in  the  wave  amplitude. 
These  singularities,  whicn  arise  along  curves  known  as 
caustics,  necessitate  the  inclusion  of  diffraction  effects, 
by  which  energy  is  allowed  to  cross  the  geometric  rays, 
leading  to  uniform  solutions  in  the  neighborhood  of  the 
singularity . 

The  combined  effect  of  refraction  and  diffraction 
was  first  studied  in  order  to  provide  locally  valid 
solutions  in  the  neighborhood  of  singularities  of  the 
geometric  optics  approximation.  Ludwig  (1966)  showed  that 
the  linear  wave  field  in  the  neighborhood  of  a  caustic  is 
represented  by  Airy  functions,  which  are  sinusoidal  in  the 
illuminated  zone  upwave  of  the  caustic,  decay  exponential  iy 
in  the  shadow  zone,  and  are  bounded  over  the  entire  domain. 
The  wave  field  in  the  vicinity  of  the  cusp  of  a  caustic  has 
been  studied  by  Smith  (1976a).  The  problem  of 
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simultaneous  refraction  and  diffraction  in  a  domain  with 
depth  inhomogeneities  was  successfully  studied  by  Berknoff 
(1972) ,  who  proposed  an  elliptic  equation  of  the  form 
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which  must  be  solved  as  a  boundary  value  problem  together 
with  appropriate  boundary  conditions.  Here,  CJ  ,  C,  and 
are  the  absolute  wave  frequency,  phase  speed  and  group 
velocity,  respectively,  and  will  be  defined  in  Chapter  2. 

A 

The  quantity  £  is  related  to  the  three  dimensional 


potential  function  ">  by 


cj>  =  P(  Cosh  k'Jn+z)  2  (/.  2.) 

Cosh  kk 

where  k  and  h  are  the  local  wave  number  and  water  deptn. 
Radder  (1979)  showed  that  (1.1)  could  be  transformed  to  a 
Helmholtz  equation  by  the  scaling 
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)  9, 


and  developed  a  parabolic  equation  approximation  for  (1.1) 
for  the  case  of  nearly  unidirectional  wave  propagation. 
Smith  and  Sprinks  (1975)  derived  (1.1)  by  means  of  Green's 
identities  applied  to  (1.2)  ,  and  found  several  higher  order 
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terms  proportional  to  the  square  and  derivative  of  the 
bottom  slope;  these  terms  will  be  discussed  in  Chapter  3. 
Booij  (1981),  using  a  method  based  on  a  variational 
principle,  derived  a  time-dependent  extension  of  (1.1) 
including  an  ambient  current;  his  equation  can  be  written  in 
the  slightly  altered  form 

,ia) 


where  £(x,y,t)  is  the  horizontal  ambient  current,  D/Dt  is  a 

A* 

total  derivation  following  the  current,  and  £  is  related  to 


i  /  ^ 

0  *  -olii  tfj'i+'z)  b 


cosh 


In  kin 


The  frequency  CT  is  the  frequency  of  oscillation  relative 
to  the  moving  domain.  Booij's  model  is  reducible  to  (1.1) 
after  neglecting  the  current  £  and  assuming  purely  harmonic 

A 

motion  (.£  =0,  where  subscript  t  denotes  tne  time 
derivative).  Booij  further  derived  a  parabolic  equation 
approximation  to  (1.4).  Equation  (1.4)  neglects  a  small 
term  which  appears  in  Booij’s  equation  due  to  an  error  in 
his  derivation. 


In  the  conclusion  of  his  dissertation,  Booi]  states 
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that  the  principle  drawbacks  to  the  model  (1.4)  include  tne 

inability  to  model  the  nonlinear  effects  of  amplitude 

dispersion  directly,  and  the  restriction  to  a  monochromatic 

wave  field.  The  second  conclusion  is  not  strictly  true, 

A 

since  allowing  ^  to  be  specified  as  a  time  dependent 
quantity  allows  for  the  specification  of  a  wave  field  with  a 
narrow  spectral  band  width,  as  long  as  a  carrier  wave  of 
frequency  CC  is  identifiable.  This  allowance  is  critical  to 
the  study  of  nonlinear  waves  in  the  context  of  the  nonlinear 
Schrodir.ger  equation,  as  discussed  below.  The  first 
conclusion  is  basic  to  the  linearization  of  the  problem. 
Booig  suggested  several  means  for  including  the  effect  of 
nonlinearity  based  on  empirical  modifications  to  the  linear 
dispersion  relation  suggested  by  Walker (1976)  and  Hedges 
(1976);  a  similar  modification  based  on  tne  correct 
nonlinear  dispersion  relation  for  Stokes  waves  will  be 
suggested  in  Chapter  3.  However,  the  present  study 
emphasises  the  inclusion  of  nonlinearity  in  the  problem 
formulation  and  model  equation  from  the  outset,  and  such  a 
model  equation  involving  a  cubic  nonlinearity  is  proposed  in 
Chapter  3. 


The  general  hyperbolic  wave  model  can  be  reduced  to 
a  first  order  hyperbolic  equation  governing  the  evolution  of 
the  amplitude  envelope  of  a  wavetrain  propagating  in  a  known 
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(or  assumed)  direction  by  substituting  the  assumed  form 


4>  ~  a6^,-0 


■f 


iU 


where  T  is  a  suitable  phase  function  for  plane  waves,  into 
the  second  order  model.  The  resulting  evolution  equation  is 
shown  in  Chapter  4  to  be  consistant  with  previous  results 
obtained  using  multiple  scale  expansions  of  the  governing 
equations.  Equations  of  this  type,  based  on  extended  forms 
of  the  Scnrodinger  equation  with  a  cubic  nonlinearity,  have 
been  shown  to  be  useful  in  modelling  effects  in  nonlinear 
wavetrains  ranging  from  predictions  of  the  stability 
characteristics  of  various  solutions  to  the  description  of 
the  evolution  of  complex  wave  patterns  in  the  presence  of 
currents  and  wind.  A  recent  and  comprehensive  review  of 
this  topic  has  been  provided  by  Yuen  and  Lake(1982).  The 
correspondence  between  the  first  order  evolution  equations 
and  the  second  order  hyperbolic  wave  equation  described  here 
points  out  the  possibility  of  modelling  complex  wave  forms 
with  narrow  spectral  bandwidth,  a  possibility  not  fully 
recognised  by  Booij. 


The  general  study  of  waves  in  a  time-dependent,  two- 
space-  dimensional  approach  involves  the  use  of  complex 
numerical  schemes  such  as  the  ADI  method.  In  this  study,  we 
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restrict  our  attention  to  steady  plane  waves.  Parabolic 
approximations  for  Stokes  waves  with  and  without  an  imposed 
current  field  are  derived  in  Chapter  4  and  used  to 
investigate  a  number  of  examples  in  Chapter  o.  In 
particular,  a  comparison  between  the  data  set  of  Berkhoff, 
Booy  and  Sadder (1982)  and  predictions  of  the  nonlinear 
parabolic  approximation  points  out  the  usefulness  of  the 
method  in  providing  estimates  of  the  wave  amplitude  in 
regions  where  ray  crossing  occurs  in  the  refraction 
approximation,  where  both  diffraction  and  nonlinearity 
become  important. 


in 

Chapter  5 

,  the  effect  of 

higher  order  t 

erms 

i nvolving 

the 

bottom 

slope  on  the 

prediction 

of 

the 

waveleng  th 

of 

edge 

waves  on  steep 

bottoms  and 

on 

the 

prediction  of  reflection  coefficients  for  waves  passing  over 
an  underwater  slope  are  investigated  using  an  elliptic 
formulation  for  steady,  linear  waves. 


CHAPTER  2.  FORMULATION  OF  THE  PROBLEM 


2.1  Governing  Equations 

In  this  study,  the  irrotational  motion  of  an 
inviscid  fluid  is  considered.  The  effects  of  surface 
tension  are  also  neglected.  The  fluid  domain  is  bounded 
below  by  a  rigid,  stationary  bottom  and  above  by  a  free 
surface,  and  is  of  infinite  lateral  extent  unless  otherwise 
specified.  A  definition  sketch  of  the  coordinate  system  and 
physical  domain  is  given  in  Figure  2.1.  The  domain  is 
allowed  to  have  an  arbitrary,  horizontally  varying 
two-dimensional  distribution  of  ambient,  quasi-steady 
current,  with  the  theoretical  restriction  that  the  ambient 
current  must  also  be  irrotational.  This  restriction  is 
seldom  satisfied  in  practice,  although  it  is  known  that  weak 
ro tational ity  does  not  effect  the  governing  equations  for 
the  geometric  optics  approximation,  as  shown  by  Mei(1982) 
and  others.  The  neglect  of  friction  also  restricts  us  to 
consideration  of  vertically  uniform  currents,  which  may  be 
considered  in  some  cases  to  represent  the  depth  average  of 
vertically  non-uniform  currents.  It  will  be  demonstrated 
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below  that  the  effect  of  friction  may  be  included  in  an  a 
posteriori  sense;  however,  it  may  not  be  conveniently 
represented  in  the  context  of  the  present  derivation. 


The  equation  of  motion  appropriate  to  the  given 
assumptions  is  Laplace's  equation  in  three  dimensions  for 
the  velocity  potential  b  : 


O 


i 

0> 


as 


o 


-  h0  —  £  ^ 


(z.  1.  i ) 


Here,  x  and  y  are  the  horizontal  coordinates,  z  is  directed 
vertically  upward  from  the  still  water  surface,  is  the 
instantaneous  free  surface  position,  and  h&  is  taken  to  be 
the  water  depth  in  the  absence  of  any  motion.  Subscripts  x, 
y,  z  and  t  will  be  taken  to  denote  differentiation  with 
respect  to  the  spatial  coordinates  and  time,  respectively. 
Kinematic  boundary  conditions  are  prescribed  for  the  bottom 
and  free  surface  and  are  given  by 


and 


v,  'nb  = 


=  - 


v.  *-• 


t  • 


>• 
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Here,  Vj,  is  3  horizontal  gradient  operator  given  by 

^  F  =  Fx  ^  +  S  *1 

where  ✓<!.,  ,  /  and  x-  are  unit  vectors  in  the  spatial 

**»  A  t 

directions.  Fluid  velocities  are  given  by 

a  *  H  <}> 

Finally,  a  dynamic  condition  is  given  on  the  free  surface  by 
the  choice 

D  C  X .  m  .  t  )  =  O  :  2 r  =  r 

i  1  ( 

where  p  is  the  fluid  pressure.  Use  of  Bernoulli's  equation 
at  the  free  surface  then  leads  to  the  condition 

<f>€  *  ±(v4>f  -o  i  5  =  7 

The  irrotational  motion  of  an  inviscid  fluid  with  uniform 
density  and  a  free  surface  and  no  applied  pressure  forces  or 
surface  tension  is  completely  described  by  (2. 1.1-4) 
together  with  appropriate  initial  and  boundary  conditions, 
which  will  be  developed  as  needed. 

The  set  of  equations  governing  the  fluid  motion  is 
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nonlinear.  in  the  present  study,  we  will  restrict  our 
attention  to  weakly  nonlinear  motions,  leading  to  the 
assumption  that  the  motion  can  be  adequately  described  by  a 
series  expansion  in  powers  of  a  small  parameter,  with 
truncation  after  a  small  number  of  terms.  However,  rather 
than  expanding  the  governing  equations  (2. 1.1-4)  directly, 
we  proceed  instead  by  using  a  variational  approach  involving 
a  Lagrangian  for  the  conservative  motion.  The  Lagrangian 
can  then  be  shown  to  lead  to  the  correct  set  of  governing 
equations  after  taking  the  arbitrary  variations  with  respect 

t 

to  tne  aependent  parameters  <J>  and  V)  . 

1  ( 

The  following  derivation  includes  the  effect  of 
interaction  between  the  wave  motion  and  the  ambient  mean 
current,  and  it  is  useful  here  to  define  several  quantities 
before  proceeding.  The  ambient  current  U(x,y,t)  will  be 
represented  ioy  tne  first  term  in  the  expansion  for  ip  ,  and 
will  be  presumed  to  be  0(1)  with  respect  to  the  small 
expansion  parameters.  Consideration  of  the  linearised 
problem  for  water  waves  in  the  presence  of  a  current  leads 
to  a  leading  order  solution  for  tne  wave  potential  of  the 
form  (see  Dean  and  Dalrymple ( 198 3) ,  sec.  3.4.5). 

.  rl 

6  s  -/i'CL  COtW  k  fh+Z.)  S.  (2-  I.  S’) 

1  S  cr 
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where 


the  phase  function. 


given  by 


The  dispersion  relation  between  the  intrinsic,  or  relative, 
frequency  J~  and  the  vector  wavenumber  k=  rf  is  given  by 


CO  =  O'  +  k '  U. 


(z.u) 


\ 

where  is  the  absolute  frequency  given  by  -  Y^  •  and 

cr  -  (ak  TAJJ»  kit  )  l 

a 


The  frequency  C"  may  be  interpreted  as  the  wave  frequency 
relative  to  a  frame  of  reference  moving  with  the  local 
current  velocity  j£. 


iiiiWmriiM—liliHillniii 
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2.2  The  Variational  Principle  and  Derivation  of  Euler 
Equat ions 

A  variational  principle  equivalent  to  the  boundary 
value  problem  given  by  (2. 1.1-4)  can  be  written  as 

sjy  L  C  4>,  <f>t  ^  ,  Vf)  ctt  cix  =0  6- 2-0 

i*, 

where  L  is  the  Lagrangian  governing  the  quantities  C  and 
)j,  t|  and  t^  are  two  arbitrarily  chosen  values  of  time,  and 
x,  given  by 

X  =  fx^l  -  X^x  + 


is  the  horizontal  position  in  the  domain.  Here  the  space 
{x,t}  may  be  considered  to  be  the  propagation  space,  and  the 

•V 

coordinate  z  may  be  considered  to  be  the  cross-space,  in  rhe 
terminology  of  Hayes(1970).  The  form  of  L  was  given  by 
Luke(1967)  as 


where  the  integration  over  z  represents  integration  over  the 
cross-  space,  leading  to  a  variational  principle  governing 
motion  in  the  coordinates  of  the  propagation  space.  The 


■tiittAliaaBiMHi 
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original  problem  (2. 1.1-4)  can  be  recovered  by  taking 

i 

variations  of  (2.2.3)  with  respect  to  an<3  P  ;  see 

1  L 

Whi tham ( 1974 )  ,  section  13.2  . 

Whitham(1967)  used  the  variational  principle  (2.2.3) 

to  derive  the  dispersion  relation  and  conservation  laws  for 

wave  action,  momentum  and  phase  for  weakly  nonlinear  Stokes 

waves  with  slow  amplitude  modulation.  The  Lagrangian  was 

treated  in  an  averaged  form  given  by 

IT 

at  (A,^)  =  L*.f  Cz;.*) 

0 

where  A  is  the  amplitude  of  the  wave  form.  Variation  of  u 
with  respect  to  A  yields  the  dispersion  relation 

&  C  <r}  k  )  +  A 1  HC k)  +  oCtH)  =  o  r) 

where 

6cVJ/?)=[;-  (u-k-iit/r*} 

is  the  dispersion  relation  (2.1.6)  for  the  linear  problem. 
Peregrine  and  Smith(1979)  have  investigated  the  effect  of 
the  nonlinear  contribution  H((T,k)  on  refraction  in  the 
vicinity  of  caustics;  this  problem  will  be  discussed  in 
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'•fi 

Chapter  6.  Conservation  laws  are  obtained  from  through 
variation  of  the  phase  function  ,  and  use  of  Noether's 
theorem . 


In  the  present  study,  the  Lagrangian  is  utilized  in 

its  unaveraged,  or  primitive,  form  in  order  to  obtain 

/ 

equations  governing  the  fast  varying  quantities  y  and 

directly;  the  resulting  equations  thus  need  no  assumption 

si, 

about  the  form  of  the  phase  function  ,  .  This  method  was 
applied  by  Booij(1981)  to  derive  a  time-dependent  mild-slope 
model  equation  for  linear  waves  in  the  presence  of  a 
current.  Booij  did  not  make  use  of  a  perturbation  expansion 
to  order  the  terms  in  his  Lagrangian;  consequently,  steady 
terms  appeared  in  the  governing  equation  for  the  linearized 
wave  motion  3nd  were  dropped  a  posteriori.  The  present 
study  provides  a  consistant  means  for  deriving  the 
linearized  equation  along  with  the  higher  order  equations  to 
be  discussed  below.  Some  use  will  be  made,  however,  of  the 
properties  of  Wh i tham ' s ( 1967 )  approach;  in  particular,  it 
will  be  argued  that  terms  which  vanish  upon  integration  over 
the  phase  do  not  contribute  to  the  governing  equation  for 
the  motion,  since  they  can  play  no  role  in  determining  the 
corresponding  conservation  laws  in  the  manner  employed  by 
Whi tham ( 1967 ) .  However,  in  order  to  strengthen  this 
argument,  it  will  be  shown  that  the  retention  of  these  terms 
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and  proper  treatment  within  the  context  of  the  perturbation 
expansion  leads  merely  to  redundant  conditions  on  the 
dependent  variables;  this  will  be  demonstrated  by  example 
rather  than  proven  conclusively. 

The  Lagrangian  L  is  manipulated  in  order  to  obtain 

1 

equations  governing  the  instantaneous  motion  of  y  and  , 
in  the  form  of  an  Euler  equation  for  (j>  and  K?  and  a 
corresponding  free  surface  boundary  condition  for  lO  .  The 

ry 

potential  y  (x,t)  is  the  portion  of  the  potential  describing 

! 

motion  in  the  propagation  space.  After  substitution  of  y 
with  known  cross  space  structure  into  L  and  integration  over 
the  depth,  L  is  of  the  form 

L*  L(<?, 


Variation  of  L  with  respect  to  Y  yields  the 
cond it  ion 


L^CSr^)  =0  •  (2.2.7) 


which  in  turn  yields  the  boundary  condition  for  Y)  after 

using  the  arbitrariness  of  £ Y)  .  Variation  of  L  with  respect 
7-  L 

to  y  then  leads  to 
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£  fr?)  +  *  LVJ  S(k$)  =  o 

(Z.Z.9 

which,  after  partial  integration,  leads  to  the  Euler 
equation 

L?  "  (LVt  ~vh-(L^  = 

where  W  may  be  interpreted  as  a  virtual  work  term  of 
arbitrary  form  related  to  slow  spatial  and  temporal 
nonperiodicity,  and  may  be  used  to  include  the  effect  of 
energy  dissipation.  This  step  was  taken  by  Booij(1981)  to 
add  a  term  to  the  mild-slope  equation  representing  the  decay 
of  wave  energy.  Dalrymple,  Kirby  and  Hwang (1983)  have 
related  Booij's  term  to  several  known  results  for  wave 
energy  dissipation  due  to  friction;  in  Chapter  6,  the 
formulation  will  be  further  extended  to  include  the  effect 
of  breaking  due  to  depth  limited  motion. 

Various  investigators  (see,  for  example,  Jimenez  and 
Whitham  (1976)  and  Ostrovskii  and  Pel inovski i ( 1972 ) )  have 
developed  methods  for  including  dissipation  effects  within 
the  variational  scheme  in  a  rigorous  manner. 


The  usual 


23 


approach  is  to  include  a  term  £  W  in  a  generalized  Hamilton's 
principle  which  may  then  be  written  as 

f  jf  L  -  ( f  wdzit  - o 

it  it 

Variations  then  lead  to  the  inclusion  of  terms  which  take  on 
the  role  of  the  term  W  in  (2.2.9).  These  approaches  have 
been  summarised  recently  by  Jeffrey  and  Kawahara ( 1 982 ) . 
Since  the  added  terms  resulting  from  the  rigorous  approach 
are  of  the  same  form  as  the  arbitrarily  added  terms  based  on 
the  allowable  nature  of  W,  we  will  not  go  further  into  these 
approaches  in  the  present  investigation. 


the 


2.3  The  Perturbation  Scheme  and  Series  Expansions 
Cor  and  ^  . 

Two  perturbation  scales  will  be  used  in 
following  derivation.  The  first  is  given  by  £  and  is  the 
usual  expansion  parameter  for  the  Stokes(1847)  wave  theory. 
By  nond imensional ising  the  governing  equations  with  respect 
to  the  wavenumber  k,  it  can  be  shown  that  £  is  equal  to  the 
wave  steepness  kA;  here,  we  will  retain  dimensional  forms 
and  employ  c  as  a  formal  parameter.  The  velocity  potential 
and  surface  displacement  are  then  given  by  an  expansion  in 
powers  of  £  corresponding  to  the  results  of  Stokes'  theory. 

A  second  formal  parameter  is  .  taken  to  be 

proportional  to  the  rates  of  spatial  modulation  of 
quantities  such  as  h,  k.  A,  etc.  This  is  equivalent  to  the 
use  of  a  multiple  scale  approach,  in  which  a  stretched 
horizontal  coordinate  X  ~  S  x  represents  the  length  scale 
over  which  variations  of  the  slowly  changing  quantities  are 
felt.  This  formalism  will  enter  the  derivation  through  the 
treatment  of  horizontal  derivatives;  for  example,  the 
gradient  of  the  depth  can  be  written  as 

^  £  /ix  ~  Vn  n  -  0{s) 

<"V  ** 

where  the  operator  is  then  an  0(1)  quantity  in  the 
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stretched  coordinate  system,  while  the  gradient  of  a 
wavelike  quantity  has  both  fast  and  slow  variations  and  can 
be  expanded  in  the  form 

Vh \\]  ~  +  "  °0)  + 

It  is  implied  that  the  operator  operates  on  the  amplitude 
of  the  functions  and  ^  ,  while  the  0(1)  derivative  is 

A 

related  co  the  gradient  of  the  phase  function  \  . 

In  the  derivation  to  be  presented  in  Chapter  3,  it 

will  be  assumed  that  the  spatial  variations  of  the  depth  h, 

wavenumber  k,  ambient  current  jJ,  and  relative  frequency  J' 

c 

may  be  characterised  by  a  variation  scale  o  of  0(d)  at  the 
largest.  This  choice  is  consistant  with  the  approach  of  Chu 
and  Mei (1970a, b)  who  neglected  the  ambient  current  U  (  and 
hence  the  relative  frequency  C  ) .  It  will  be  shown  that  a 

C“  - 

shift  of  scaling,  with  o  taken  to  be  of  0(c  ') ,  then 
reproduces  the  mild  slope  formulation  of  such  authors  as 
Ber khof f ( 1972 )  and  Booij(1981). 

In  contrast  to  the  relatively  fast  spatial 

variations  allowed  by  the  choice  of  4  ~  ,  we  will  restrict 

time  variations  of  the  quantities  h,  k,  U,  and  'C  to  be 

rS 

proportional  to  0 (  o  ) }  this  will  have  the  effect  of 
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excluding  time  variations  of  the  physical  domain  from  the 

derivations.  This  choice  is  consistant  with  an  assumption 

that  the  time  dependence  of  the  physical  domain  is  related 

to  tidal  variations  at  the  fastest,  and  corresponds  (  in  a 

multiple  scale  approach)  to  the  choice  of  a  stretched 

■» 

coordinate  t  =  t  for  variations  other  than  those  of  the 

principle  wave.  For  example,  if  the  incident  wave  is 
characterised  by  a  period  of  10  seconds  and  a  steepness 
c=0.1,  then  the  slow  variations  in  t%  will  occur  over  a 
time  scale  of  0(  20  hours  ) . 

In  the  perturbation  expansion  schemes  standardly 
used  for  deriving  nonlinear  evolution  equations  (  see  Davey 
and  Stewar tson ( 1974 ) ,  Chu  and  Mei (1970b),  and  Yue(1980)  for 
typical  derivations  related  to  Stokes  waves  in  intermediate 
water  depth  ),  the  modulation  parameter  £  is  typically 
related  directly  to  the  Stokes  expansion  parameter  6  ,  with 

disparities  in  scale  maintained  by  the  use  of  multiple 
scales  in  the  propagation  space.  Then,  the  expansions 
proceed  in  powers  of  the  single  parameter  i  .  Here,  the 
parameter  £  will  not  be  related  in  size  to  the  parameter  L 
until  the  point  is  reached  where  specific  results  for 
certain  choices  of  scale  relations  are  derived.  This 
approach  has  certain  specific  advantages  over  the  standard 
single  parameter  approach,  in  that  it  does  not  become 
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necessary  to  Four i er-expand  the  boundary  value  problems  at 
each  order,  since  this  need  is  a  direct  result  of  having 
lumped  the  expansion  parameters  in  the  first  place.  Thus, 
terms  which  are  implicitly  related  to  higher  order  in  £  at 
a  given  order  of  6  do  not  become  explicitly  tied  to  terms 
at  higher  order  in  £  . 

In  order  to  obtain  the  desired  results  for  weakly 

nonlinear  waves,  it  is  necessary  to  carry  the  calculation  of 

y  t 

the  Lagrangian  L  to  0(£  ),  based  on  expansions  of  and  >'< 

to  0(6)*  We  will  impose  the  restriction  that  the  smallest 

terms  at  each  stage  of  the  expansion  have  exponents  of  i 

and  £  totalling  four;  thus  terms  in  L  of  0  (  s'1)  must  retain 
r  i  y 

terms  to  0 (  b  ),  while  terms  of  0(6  )  retain  no  modulation 
terms.  This  will  have  the  effect  of  excluding  slow 
modulations,  and  thus  any  consideration  of  scale,  from  the 
nonlinear  terms  found  below.  This  restriction  is  consistant 
with  che  theories  of  Chu  and  Mei (1970b)  and  subsequent 
investigators,  with  the  exception  of  Dys the ( 1 97 9 )  ,  who 
includes  the  lowest  order  modulation  effects  in  the 
nonlinear  contribution. 

One  drawback  of  the  Lagrangian  approach  is  that  it 
provides  no  information  about  the  correct  form  of  the 

i 

expansions  of  <p  and  .j  to  be  used,  both  in  terms  of 
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consistency  with  the  ordering  parameters  and  in  terms  of  the 
structure  of  the  solution  in  the  cross-space;  this  has  led 
to  several  misunderstandings  in  the  past.  Wh i tham ( 1 96 7 ) 
assumed  a  one-dimensional  (in  x)  expansion  of  (j?  and 
given  by 

=  6  fio  4  +  £  XJio  4  *  y  "  '^-3.!  a) 


0 


where 


x 

Tib 


cask  k  Ch+z) 

I  {_  1 

CoJ77 


css*  h  +:£S) 
S/AJf}  '  krt 


and  where  £>  ,  T  ,  and  b  are  0(c).  The  resulting 

/ 

conservation  equations  were  subsequently  found  to  be 
inconsistant  with  multiple  scale  expansions  based  on 

r> 

modulations  of  0(a  )  0 ( c  ),  as  shown  first  by  Chu  and 
Mei (1970b) ,  and  it  was  surmised  that  Whitham's  theory  was 
basically  incomplete.  Wh i tham ( 1 97 0 )  showed,  by  means  of  a 
multiple  scale  expansion,  that  his  results  were  consistant 
with  the  modulation  rate  ?  £  .  Finally,  Yuen  and 
Lake(1975),  for  the  case  of  deep  water,  showed  that  the 
inclusion  in  Whitham's  approach  of  the  proper  forms  of  ^ 
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i>  =  b‘  C<  t)  t  %z&h 

k  O  i  o  ^  /-  >  *  - 


// 


Cz.S-^, 


where 


cp 

TO 


/  i  ,  V  i  » 

i  ”» *?/  y  ?« 


-  (.  !n0  +  i)  V,  /l„  -y, 


/  -»  •?  //Li 

u.  ■»■  v°y 


The  quantities  f j ^  and  > .  are  given  by 


-£  =  j?  c  cos'n  k  i  n 

cosk  kh 


r  u  /,  ,  "s  j_2  ^  '  '  /  '  ■'  /  I  /  / ,  _ 

-r .  =  g,  /i— t;  J.x/7  -  Rsi  -T^jun  Rn  coshR^Zy 

J  !<■ 


L  kh 


Cos  J>  kh 


l  =f  kZi'n  +  z)  -  Celt)1?  ct>i'*kC»+z) 

J  1  j  cesUMh 


C"  s.r; 


and 


*■  i  ~  j  <*' 

* 


)  *** 


£-Vh/?  /“.I.o" 


ifc* 


where  h  is  the  0(1)  mean  water  level  given  by 
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-  f, 


0  & 


(2,3.7) 


The  forms  of  the  wavelike  components  of  Q  and  O  are 

1  L 

essentially  similar  to  those  given  by  Chu  and  Mei(1970a), 
with  the  exception  that  the  absolute  frequency  is 

replaced  everywhere  by  the  relative  frequency  (J  . 


Other  terms  included  in  (2.3.3)  are  b0  and  b ^  , 

which  are  slowly  varying  surface  terms  related  to  the 

ambient  mean  water  level  and  the  wave- induced  surface 

1  /  2. 

fluctuation;  p  ,  which  is  the  potential  for  the  0(6  ) 

v'  V 

wave-induced  current,  and  vQ  and  0^  ,  which  are  Bernoulli 
constants  for  the  mean  flow  at  each  order. 


In  order  to  place  the  present  derivation  within  the 
context  of  the  implied  or  explicit  scaling  assumptions  of 
various  studies,  a  list  of  the  various  choices  is  given  in 
Table  2.1. 


Highest  order  Modulation  I  Modulation  I  Largest 
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CHAPTER  3.  THE  LAGRANG IAN  AND  EQUATIONS  GOVERNING 


WAVE  MOTION 


3.1  Introduction 

The  Lagrangian  L  governing  the  motion  of  waves  with 
imposed  current  field  JJ(x,t)  is  constructed  by  substituting 
tne  perturbation  expansions  (in  the  small  parameter  &  )  for 
and  t  given  by  (2.3.3),  into  (2.2.3)  and  performing  the 
integration  over  depth.  The  Lagrangian  is  then  expanded  as 
a  series  in  the  small  parameter  &,  yielding  the  expression 

L=LC  +  6  L(  +  ^  +  *3L5  +  6\  '3.  3.:) 

Then,  consistent  with  the  restriction  that  the  exponents  of 
c  and  £  total  no  more  than  four,  we  may  restrict  L*  to 
0(  i  )  and  L  ^  to  0(1).  L ^  will  not  contribute  any  new 

information,  and  may  be  truncated  to  0(1)  in  S  for 
simplicity.  Equations  governing  the  various  terms  in  <j>  may 
then  be  obtained  by  varying  the  term  by  the  member 

of  the  perturbation  series  for  .  Since  the  perturbation 
method  is  somewhat  unique,  an  example  of  its  use  in  a 
simpler  problem  is  shown  in  Appendix  B. 
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3.2  The  Lagrangian  for  Stokes  Waves 


The  assumed  forms  for  •{>  and  l ^  are  given  by 
(2.3.3),  based  on  the  derivation  for  the  large  current  case 
giver,  in  Appendix  A.  Individual  contributions  to  the 
Lagrangian  are  given  by  gz,  and  V|  ;  the  subsequent 
calculations  are  lengthy  and  are  outlined  in  Appendix  C.l. 
The  results  may  be  written  in  the  expanded  form  (3.1.1);  the 
individual  contributions  are  given  by 


t-o  -  f  a  f 4  ' X> 4' )  ]  -1'  t 

T-  WX  LX.  )  4? 


') 


where  h=n0+  b0  and  wnere  <4  and  are  t0  oe  expanded  in 

the  forms 


<k 


.  / 

<?0  *• 


r 


4* 


u 


and 


e> 


u 


and  we  assume  that  the  time  variation  of 


l  '  p* 

r t  is  0(S); 
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•  / 

where  xi  is  also  an  integral  expression  evaluated  in 
Appendix  D; 


and 


3.3  The  Governing  Equations 


We  proceed  by  determining  the  Euler  equations  for 
each  order  in  <£  ,  following  the  method  outlined  in  Appendix 
B.  First,  by  varying  L0  with  respect  to  b0  ,  we  obtain  the 
expression 

<j£  .  <j>„'  -  x  »  o.  0(T)  (3  5./) 


By  requiring  that  b0  vanish  in  the  aosence  of  a 

V 

set  =0.  Then,  varying  L0  with  respect 

performing  the  partial  integration  over  the 
space,  we  obtain  the  continuity  equation 


current,  we 

i 

to  and 

propagation 


1  =  0-  0(V)  a.3.3) 


By  taking  the  horizontal  gradient 
ODtain  a  momentum  equation  for  the  0(1)  mean 

*&■%)'£  *  *o 


of  (3.3.1),  we 
flow: 

^5  3.3) 


where  j,J  -  •  Recalling  from  (3.2.1)  tnat  b&  and  are  to 

r  }  r* 

be  expanded  in  orders  of  a,  and  tnat  ■^■'w0(  »  )  for  the  zeroth 
order  flow,  we  find  that 


-  OfS1) 
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ybj  =  -{('Z?.')1  +  O  (s1)  (itr) 

to  leading  order.  All  contributions  to  Lj  and  L1  involving 

,  / 

T/j,  ( h '/,,£)  may  therefore  be  dropped.  The  contribution  of 
•  '/ 

in  L,  is  thus  reduced  to  a  small  correction  to  the 

operator  D/Dt  and  the  0(1)  mean  water  level  h,  and  enters 

,  / 

the  overall  equation  in  tandem  with  .  It  should  be  noted 
that  in  ,  which  is  expanded  only  to  0 (  £  ),  the 
contribution  from  b6  and  <pa  in  h  and  D'/Dt  should  be 
dropped;  however,  retaining  them  introduces  only  a  small 
error.  It  is  apparent  from  (3. 3.4-5)  that  the  current  is 
effectively  steady  to  leading  order,  and  the  smallness  of  o0 
justifies  the  neglect  of  time  derivatives  of  tne  f  functions 
in  i  . 

i 

Moving  to  0(6  ),  we  note  that  the  average  of  L,  over 

the  phase  for  any  choice  of  temporally  periodic  function  for 

i p  vanishes;  i.e., 

T  Zt 

o 

Thus  L(  does  not  contribute  in  any  significant  way  to  che 
governing  and  conservation  equations,  following  the  results 
of  Whi tnam ( 1967) .  As  a  test  of  this  hypothesis,  we  need  to 
show  that  variations  of  L  t  produce  only  redundant 
conditions  on  be  and  ^  .  At  0(  6  )»  our  lowest  order 
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r~i  t  / 

unknowns  are  ^  and  <^>(  ,  since  we  can  assume  that  bQ  and  ^ 
are  specified  by  (3. 3. 1-2).  Varying  L;  with  respect  co  tj, 

/y 

gives  (3.3.1)  again.  Varying  L(  with  respect  to  Z't  after 
dropping  the  0(  S^)  contribution  of  H'  and  J1  leads  to  the 
condition 


[&c  -  S-  £ 

which  is  accurate  due  to  our  assumptions 
scale  of  variation  of  the  current. 


about  the  time 


Moving 

respect  to  the 
fj^i  variation 
condition: 


to  0  (  )  ,  we 
still  unknown 
leads  to  the 


take  variations  of  with 

/T 

quantities  and  y,  .  The 

well  known  surface  boundary 


(3.3.7) 


where  ' =  1.  The  variation  with  respect  to  and  the 

subsequent  partial  integration  are  tedious  and  are  presented 
in  Appendix  C.2.  The  resuiting  second  order  wave  equation 

/y 

for  ^  is  given  by 


oe  ^  '  Dt 


-  ^  ^  J 


(3.3.?) 


where  F  is  a  complicated  coefficient  given  below  and  derived 
in  Appendix  C.2.  The  term 

is  identical  to  the  higher  order  terms  found  by  Smith  and 
Sprinks  (1975);  however  it  is  not  the  only  correction  at 
0  (  <5  )  involving  the  bottom  slope.  The  term  F  is  given,  for 
the  general  case  of  currents  and  depth  varying  at  0(S  ),  by 


s,  =  kvh  (Vm) 

I  >  jCZ. - 2  I  ( - 

'  k'A 


=  VL  i  fe  T  ^7  j 


.  U  i 

ft  r 


h  *  in  .Vs,)/:** 


3.3. 


are  0(S  5  quantities,  and 
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are  0($  )  quantities  appearing  in  products.  The  quantities 

v/  *✓  *"•* 

J|  -  c'^  and  <?'  arise  from  integrals  of  the  normal  mode 
structure;  explicit  values  are  given  in  Appendix  D.  The 
term  can  be  expanded  explicitly  as 

o(  =  +  +  $"3^4  k  ~^h\  ^  T‘  (^3.3,/3) 

and  £,  -  bf-  are  also  given  in  Appendix  D. 

Equation  (3.3.3)  is  an  equation  for  linear  waves  in 
a  domain  with  depth  and  currents  which  vary  at  0(S  );  as 
such,  it  serves  as  the  first  higher  order  correction  to 
Booi j ’ s (1981)  model,  which  is  essentially  similar  to  (3.3.8) 
after  neglecting  F  and  . 

Moving  to  L  3  ,  variation  with  respect  to  the  unknown 

second  order  free  surface  (b  ^  ^ z  )  yields  (3.3.7)  as  a 

/■*** 

redundant  ondition.  Varying  L3  with  respect  to  ^  yields  a 
second  governing  equation  similar  to  (3.3.8);  this  equation 
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can  be  shown  to  be  satisfied  at  lowest  order  by  waves 
governed  by  the  linear  dispersion  relationsnip  (2.1.6),  and 
is  therefore  also  a  redundant  condition. 


Variation  of  L ^  with  respect  to  (•^l+b1)  yields  the 
free  surface  condition, 


+  <p6  •  Vh  +  1  Ga  rj 

-  H  ^  //  A*-  ^  -5  11  /N<  t 

4  i!  24  +  ^  fa  +  4a  ?,  =0 

Dt  z 


Dt 


(l  2.  a) 


which  specifies  the  complete  second  order  free  surface.  The 
mean  surface  b^  can  be  extracted  from  this  expression  by 
taking  the  mean  over  a  wave  period.  After  substituting 
(3.3.7),  this  operation  results  in 


aOi  +  ‘w  Oj.  - 

o  1  Ot 


///,  , 
■S  *  Dt  > 


z. 


i^fa 


*3 


/'  O'  1 

5ft,  o,  (3. 3.  iff) 

■*» 


Again,  by  requiring  that  b^  vanish  in  the  absence  of  waves, 

Y  I  ' 

we  set  0.  Then,  varying  with  respect  to  $ l  leads  to 

the  expression 


-  ?J  *  o 


(jw) 


Jd 
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Taking  the  mean  over  a  wave  period,  (3.3.16)  reduces  to 

This  is  the  continuity  equation  for  the  second  order 

•  / 

wave-induced  current  vi,^.  Writing  the  R.H.S.  of  (3.3.15)  as 

2 

and  eliminating  b^_  from  (3.3.15-17)  leads  to  a  forced 
wave  equation  for  the  second  order  mean  flow,  given  by 


This  equation  is  equivalent  in  spirit  to  the  forced  wave 
equation  given  by  Davey  and  Stewartson  (1974)  after 
including  the  effect  of  0(1)  currents;  the  connection  will 
be  seen  more  clearly  after  substituting  for  !3  and  7J,"  (  "A  4j ) 
based  on  plane  wave  approximations. 


Variation  of  L  ^  with  respect  to  G>t  leads  to  the 


governing  equation 


•Vf{AV.Ji  -7*-UVk^'J  5.? 


^  4  -© 


(?.  3.  n; 


(3.3.16)  and  (3.3.19)  may  be  combined  by  eliminating 


( /^t+bL)  ;  the  resulting  equation  governs  the  second  order 
potential  .  Equations  (3.3.2),  (3. 3. 7-8),  (3.3.18)  and 
(3.3.19)  represent  the  complete  set  of  governing  equations 
for  the  Stokes  wave  train  to  0{6L),  with  modulation  rates  of 
0(^6)  allowed  for  all  quantities  with  the  exception  of 
time-variations  of  . 


3.4  The  Mild  Slope  Approximation  and  Some  Plane  Nave 


Results 


Reduction  of  the  previously  derived  models  to  the 
case  of  very  slowly  varying  (0(b)  )  currents  and  depth  can 
be  obtained  by  truncating  the  equations  to  0(S  ).  This 

.  #  i 

removes  the  effects  of  the  oi  ■  and  o  from  d>  ,  and  leads  to 

jl  • a  > 

the  neglect  of  the  coefficients  F  and  h*'  in  (3.3.8)  .  (  It 

is  noted  that,  in  some  instances,  it  may  be  desirable  to 
retain  the  fast  modulation  rate  for  the  wave  amplitude, 
contained  in  the  coefficients  |$3  and  ,  due  to  the 

influence  of  nonlinear  effects.)  The  first  order  wave 
equation  becomes,  after  substituting  for  the  coefficients 


from 

Appendix  D 

D £ 

i)$f  Vki  !  +tr'C 

ii 

"5* 

i 

(2,V./) 

which 

i  is  equivalent  to  Booij's 

model  after 

neglecti 

ng  an 

0(  S") 

term  which  appears  due  to 

an  error  in 

Booi  j  '  s 

free 

surface  boundary  condition. 

The  terms  at 

0(  c’) 

are 

unaltered,  since  modulations  we 

re  not  consid 

ered  at 

that 

order 

•  • 

Several  expressions  from 

section  (3.3) 

may  be 

cast 

in 

more  familiar  form  by 

introducing 

plane 

wave 

approximations  for  the  terms  and  <?(  ,  obtained  from 

1  “  *  ’  1  1  * 

Appendix  A.  The  term  Vh  •  j  ^  7h  f(  ]  may  be  reduced  to  the  usual 
expression  for  the  wave-induced  mass  flux 

where 

£-i<V 

In  the  absence  of  a  mean  current  U,  (3.3.17)  becomes  the 
familiar  continuity  equation 

♦  V  {<>.?»  *  *§k] -o  »*»> 

found  originally  by  Whitham (1962)  and  Longuet-Higg ins  and 
Stewart ( 1962) .  The  extension  to  the  case  including  a 
current  is  given  by 


C3.  i.i) 


which  reflects  the  interaction  between  the  mass  flux  of  the 
0(1)  current  and  the  small  correction  to  the  total  depth 
given  by  b^ . 


In  (3.3.18),  the  quantity  may  be  estimated  by  the 


plane  wave  approximation 


In  a  quasi-steady  wave  field  with  small  spatial  variations, 
bz  is  therefore  given  by 

b  *  bt  ’  ~  1 

which  is  recognisable  as  the  steady  0(  d1)  wave-induced 
setdown,  given  by  Longuet-H igg ins  and  Stewart ( 1964) .  The 
extension  to  the  unsteady  case  with  current  is  given  by 


where  b1$  is  given  by  (3.4.6).  The  forced  wave  equation 
(3.3.18)  may  then  be  simplified  by  the  use  of  (3.4.2)  and 
(3.4.5)  to  give 


A  form  of  this  equation  applicable  to  water  of  constant 
depth  in  the  absence  of  an  0(1)  current  is  given  by 
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given  by  Yue(1980)  and  in  slightly  different  form  by  Davey 
and  Stewartson ( 1974) . 


* 


3.5  The  Weakly-nonlinear  Equation  for  <£j 

In  the  preceeding  sections,  a  method  has  been 
described  by  which  linear  equations  governing  the  evolution 
of  the  first  two  terms  in  the  Stokes  expansion  for  (j>  ,  given 
by  £  and  ,  are  obtained.  This  metnod  could  be  carried 
to  higher  orders  after  allowing  perturbation  of  the 

3 

wavenumber  k  to  account  for  secular  terms  arising  at  0(6); 

however,  the  algebra  involved  becomes  exceedingly  tedious 

for  the  general  theory  being  developed  here.  The  results 

described  so  far  allow  for  the  calculation  of  the  0 ( t  ) 

linear  wave  field,  with  subsequent  calculation  of  the  0(ei) 

forced  harmonic  following  based  on  the  result  for  the  linear 

fundamental.  Truncation  at  this  level  provides  results 

consistant  with  the  Stokes  solution  to  second  order,  with 

the  waveform  characterised  by  short  crests  and  broad  troughs 

due  to  the  second  order  correction.  However,  the  results  do 

not  include  the  effect  of  amplitude  dispersion  on  the  phase 

2 

velocity;  this  result  would  arise  at  0(  6  )  through  the 
choice  of  the  first  perturbation  correction  to  the 
wavenumber  k. 

One  method  of  correcting  this  deficiency  could 
proceed  by  calculating  k  based  on  a  nonlinear  form  of  the 
dispersion  relation  (2.1.6).  Booij(1981)  suggested  two 


empirical  modifications  to  the  linear  dispersion  relation 
based  on  suggestions  by  Walker  (1976)  and  Hedges ( 1976 ) ; 
however,  it  is  simple  enough  to  use  the  dispersion  relation 
from  nonlinear  theory  directly,  thus  dispensing  with  the 
empirical  approximations.  By  averaging  a  Lagrangian  similar 
to  (3. 3. 1-5),  neglecting  the  0(Ji)  terms,  Whitham  obtained  a 
dispersion  relation  given  by  (including  the  current  U) 

(to-  (&+Vk  &)■'&)  =  qkTAvbkh  (l+tk!MfD)  (3.  S’.  l) 

where 


D  « 


1 

COSrt 


hkn  —  iTAuin  Rh 
ft  Hun  H  kh 


i3  S. 


The  relation  (3.5.1)  then  provides  a  means  by  which  the 
value  of  k  may  be  determined  based  on  local  values  of  the 
amplitude  A.  Inclusion  of  this  effect  in  the  linear  model 
(3.3.8)  requires  an  iteration  in  order  to  obtain  the  initial 
guess  for  the  updated  local  value  of  A,  after  which  the 
corrected  value  of  k  is  determined.  This  method  was 
suggested  by  Booij  as  a  means  for  incorporating  the 
empirical  nonlinear  dispersion  relations  of  Hedges(1976)  or 
Walker(1976)  in  the  linear  model;  and  will  be  employed  in 
section  6.9  where  Stokes  theory  is  not  appropriate  to  the 
description  of  the  problem. 


In  the  recent  literature  on  nonlinear  water  waves, 
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the  most  successful  approach  to  including  the  effect  of  low 
order  nonlinearity  in  the  basic  governing  equation,  and  to 
producing  tractable  computational  models,  has  been  to 
include  the  nonlinear  effects  in  the  form  of  nonlinear  terms 
added  to  the  linear  governing  equation.  The  added  terms, 
proportional  to  the  third  power  of  the  local  amplitude, 
account  for  the  effect  of  amplitude  dispersion  provided  by 
the  nonlinear  dispersion  relation  (3.5.1)  or  its 

fast-modulation  counterpart. 

We  proceed  by  deriving  the  third  order  correction 
terms  to  be  included  in  the  linear  governing  equation.  This 
development  leads  to  the  correct  inclusion  of  amplitude 
dispersion  in  the  linear  wave  equation,  which  was  only 
suggested  on  approximate  grounds  by  Booij (1931) .  in 
addition  the  second-order  hyperbolic  form  of  the  governing 
equation  allows  for  the  study  of  reflections  in  an  arbitrary 
domain,  which  has  been  investigated  only  occasionally  in  the 
context  of  the  nonlinear  Schrodinger  equation  approach  (  see 
for  example  Smith  (1976b)  and  Newell  (1978)  ). 

It  is  sufficient  to  consider  a  reduced  form  of  L 

given  by 

L  -  Co  C S'})  -  t-u  (Z.S.Z) 


Here,  we  have  dropped  terms  of  0  ( £  *)  from  for  the  sake  of 
computational  simplicity;  these  terms  are  unaltered  by  the 
derivation  and  may  be  replaced  subsequently.  Also,  we  drop 
Lj  entirely  on  the  grounds  that  its  average  over  the  phase 
is  zero,  implying  that  it  has  no  effect  on  the  conservation 
equations;  alternatively,  we  can  argue  that  variation  of  L  ^ 

A# 

with  respect  to  <(r  in  the  following  derivation  leads  to 

-  2  K 

terms  which  are  steady  or  proportional  to  e  ,  and  are 
asynchronous  with  the  linear  wave  component;  they  therefore 
do  not  contribute  to  the  motion  of  the  fundamental  wave. 
However,  ,  when  varied  with  respect  to  $t  and  , 

•  Z  *  7  y  /  - 1  ^ 

produces  terms  proportional  to  e  and  e  ;  the  first 

set  of  terms  is  then  retained  as  the  principle  nonlinear 

r~> 

correction  to  the  governing  equation  for  ot  . 

First,  the  reduced  Lagrangian  (3.5.3)  is  varied  with 
respect  to  <y,  ,  and  the  resulting  expression  is  partially 

integrated  to  give  the  nonlinear  counterpart  of  (C.2.1): 


]  -c -s  v'3 .r. 
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Varying  L  with  respect  to  and  dividing  by  g  leads  to  the 
nonlinear  free  surface  boundary  condition: 


fl 
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In  the  equations  (3. 5.4-5),  bt  is  given  by  (3.4.7). 
Differentiating  (3.5.5)  and  substituting  in  (3.5.4),  we 
obtain  a  nonlinear  governing  equation  for  C>  ,  with 
nonlinear  terms  including  the  first  and  second  order 
fluctuating  surfaces  ^  and  : 
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Here,  nonlinear  terms  involving  (vrt'U)  have  been  dropped  due 
to  the  assumptions  about  modulation  scales.  Equation 
(3.5.6)  may  be  written  entirely  in  terms  of  the  potential  J 
by  the  laborious  substitution  of  the  governing  equations  for 

/V 

.,'i  ,  /;t ,  and  ^  r  use  of  approximate  forms  of  these  equations 
based  on  the  local  neglect  of  modulations  does  not  alter  the 
level  of  accuracy  of  the  nonlinear  terms  in  (3.5.7).  We  will 
proceed,  instead,  to  a  much  reduced  form  of  the  nonlinear 
equation  by  substituting  the  plane  wave  results  of  Appendix 
A  in  the  nonlinear  terms  of  (3.5.6).  This  approximation 

7 

would  lead  to  errors  of  0(£  )  in  the  direct  calculation  of 
nonlinear  standing  waves,  where  the  amplitude  modulation 
cannot  be  considered  to  be  small;  however,  many  of  the 
examples  below  are  for  quasi- unidirectional  wave 
propagation,  and  the  substition  made  here  is  entirely 
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acceptable.  We  also  neglect  derivatives  of  the  amplitude  A 
during  the  substitution,  which  would  lead  to  the  inclusion 
of  inconsistantly  small  terms  of  0(£3,  $  )  .  Substituting 
(3.4.7)  and  the  plane  wave  approximations  into  (3.5.6)  leads 
to  the  equation 
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and  D  is  defined  in  (3.5.2).  Substituting  for  the  integrals 
Ggg  and  ,  dropping  the  asynchronous  terms  proportional 
to  e  and  rearranging,  we  arrive  at 
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This  represents  the  mild-slope  equation  for  ,  in 

the  presence  of  a  slowly  varying  ambient  current  U ,  correct 

•*  / 

to  0(  6*)  .  The  added  terms,  proportional  to  |A,  '  and  , 
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represent  the  direct  effect  of  amplitude  dispersion  and  the 

/ 

interaction  with  the  wave  induced  current  Vj,  y'i  > 
respectively.  Note  that  must  also  be  solved  for  in  the 

general  case,  using  the  results  of  section  3.4. 

The  fast  modulation  corrections  to  the  third  order 
model  equation  (3.5.9  )  may  be  obtained  by  including  the 
0  (  i  ,  >*)  terms  from  (3.3.8),  to  give 
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Equation  (3.5.10)  is  a  principle  result  of  this 
investigation . 


In  order  to  deliniate  the  connection  between 
(3.5.10)  and  the  results  of  previous  investigators,  we 
remark  that 


C01  *  U>*  +  co* 
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where  and  are  given  by  Chu  and  Me.*.  (1970),  and 

correspond  to  the  second  order  frequency  corrections  due  to 

amplitude  dispersion  and  change  in  mean  water  level, 

,  j 

respectively.  In  addition,  -fi.  equivalent  to  Chu  and 

•  • 

Mei ’ s  while  the  remainder  of  (3.5.10)  corresponds  to 

y  *T* 

their  CJZ"  and  :-z  ,  after  accounting  for  the  inclusion  of  the 
fast  wavelike  motion. 


CHAPTER  4.  EQUATIONS  GOVERNING  THE  AMPLITUDE  A  IN 
QUASI  ONE-DIMENSIONAL  PROPAGATION 

4.1  Introduction 

The  second  order  hyperbolic  models  (any  of  (3.3.8), 

(3.4.1),  (3.5.9)  or  (3.5.10)  )  represent  a  very  general 

approach  to  che  solution  of  wave  propagation  in  a  domain 

with  correctly  posed  initial  and  boundary  conditions; 

however,  the  existing  numerical  solution  methods  for  this 

type  of  problem  in  an  arbitrary  domain  (in  particular, 

finite  element  or  finite  difference  methods)  are  costly  and, 

in  the  case  of  finite  element  programs,  require  a  nontrivial 

programming  effort.  in  several  applications,  one  or  both  of 

two  simplifying  assumptions  may  be  appropriately  utilized. 

First,  the  wave  amplitude  may  be  assumed  to  be  steady  in 

time;  thus  the  only  contribution  of  terms  like  ^  Ti  /  }  D  is 

/"*■» 

the  0(1)  phase  contribution  This  approximation  may 

be  used  in  most  applications  of  the  linear  wave  equations 
unless  it  is  desirable  to  study  time  variations  of  the  wave 
field  due  to  unsteady  incident  wave  amplitude,  representing 
the  effect  of  wave  groups  or  transient  disturbances.  For 
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nonlinear  applications,  use  of  this  simplifying  assumption 
implies  the  tacit  assumption  that  nonlinear  interactions 
between  waves  in  a  multicomponent  system  will  not  lead  to 
the  onset  of  significant  amplitude  modulations,  such  as 
Benjamin-Feir  instabilities,  over  the  length  and  time  scales 
being  modelled;  this  assumption  is  problematic  and  has  not 
been  fully  investigated  for  the  case  of  an  initially 
unmodulated  plane  wave.  Equations  based  on  the  use  of  this 
assumption  in  order  to  reduce  the  second  order  hyperbolic 
system  to  an  elliptic  system  will  be  given  in  section  4.2. 

An  independent  assumption  of  quasi-one-dimensional 
wave  propagation  can  be  made  a  priori,  with  the  resulting 
restriction  that  the  reflection  of  waves  counter  to  the 
principle  direction  of  propagation  cannot  be  considered  in 
the  same  computational  step.  In  this  case,  the  choice  of  x 
as  the  principle  propagation  direction,  with  the  concurrent 
choice 

S'  ( fe>°! 

leads  to  a  general  set  of  equations  in  {x,t}  where  the 
dependence  in  y  is  related  to  oblique  modulations  of  the 
amplitude.  In  cases  where  refraction,  and  resulting  bending 
of  wave  rays,  is  important,  the  y  dependence  is  used  to  take 
up  errors  resulting  from  k  not  being  truly  alligned  with  x. 
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Time  dependent  equations  based  on  this  approximation  are 
given  in  section  4.3,  and  we  study  the  connection  between 
the  general  fast  modulation  model  studied  here  and 
previously  derived  models  based  on  the  nonlinear  Schrodinger 
equation,  which  are  related  to  subsets  of  the  general 
problem . 

Finally,  in  section  4.4,  an  approximation  method 
analogous  to  the  splitting  matrix  approach  {see,  for 
example,  Corones  (1975),  McDaniel ( 1975) ,  and  Radder(1979)  ) 
is  used  to  obtain  the  parabolic  approximation  to  steady, 
quasi-one-dimensional  propagation,  based  on  the  elliptic 
models  of  section  4.2.  We  also  demonstrate  the  connection 
between  the  parabolic  equations  and  the  scaling  assumptions 
neccessary  for  the  recovery  of  these  approximations  from  tne 
initially  time  dependent  models  of  section  4.3.  The 
coupling  between  forward  and  back-scattered  wave  components 
will  be  neglected. 
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4.2  Time-independent  Elliptic  Models 

Before  proceeding  with  the  reduction  of  the 
hyperbolic  model  to  an  elliptic  form,  we  choose  a  form  for 
the  wave  energy  dissipation  function  W, : 


(H.  2.l) 


which  is  slightly  different  than  the  form  assumed  by 
Booij(1981)  and  reduces  to  the  more  carefully  derived  form 
found  by  Jimenez  and  whitham ( 1976) ,  among  others.  A  method 
for  relating  w  to  several  known  dissipation  mechanisms  is 
given  in  Appendix  G.  Substituting  (4.2.1)  into  the  full 
nonlinear  model  (3.5.10)  and  making  the  assumption 
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reduces  (3.5.10)  to  the  elliptic  form 
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where  we  have  used  the  fact  that  =0  in  3  stea<^y  wave 
field.  The  corresponding  mild-slope  approximation  for 
linear  waves  is  given  by 


-2 +  u. (Ur-%  ?, )  -4  Ya -V.  ?, ) ' K •  C 5 ) 
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which  is  similar  to  equation  (3.23)  in  Booij(1981). 
Neglecting  the  0(1)  current  u  reduces  (4.2.4)  to  the  form 
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which  is  the  mild  slope  model  of  Ber khof f ( 197 2 , 19 76 ) 
extended  to  include  wave  energy  dissipation. 


The  time  independence  applied  to  the  second  order 
results  of  section  3.4  leads  to 


and 


Vu 
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where  (4.2.6)  is  a  statement  that  the  total  second  order 
component  of  the  mean  volume  flux  is  nond i vergent . 
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Equations  (4.2.4)  and  (4.2.6)  must  be  solved  as 
pair  in  the  general  case. 


a  coupled 
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4.3  Time-Dependent  Equations  for  A 


In  this  section,  we  use  the  assumed  form  for  Cjj  : 


4>, 
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to  derive  a  time  dependent  governing  equation  for  A.  The 
resulting  equation  for  one-directional  wave  propagation  is 
in  the  form  of  a  modified  cubic  nonlinear  Schrodinger 
equation,  and  subsets  of  the  general  equation  can  be 
directly  related  to  the  governing  equations  arising  in 
previous  studies  of  the  dynamics  of  nonlinear  waves. 


The  most  general  form  of  the  governing  equation  is 
given  by  (3.5.10),  and  is  rewritten  here  for  convenience. 
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The  term  (si  +F)  is  given  in  section  (3.3).  The  equation  also 
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contains  a  dissipation  term.  In  addition  to  the  wave 
equation,  we  will  also  need  the  lowest  order  energy 
equation,  derived  from  the  results  of  Appendix  E: 


At  *  (C,~u)-VhA  *  jf) 


A  'O 

(y.i.i) 


The  reduction  to  the  cubic  nonlinear  Schrodinger  equation 

( NLS )  is  made  by  substituting  (4.3.1)  into  (4.3.2)  and 

i.'f' 

cancelling  the  factor  eA'  .  The  0(1)  current  0  may  be 
written  in  its  x  and  y  components  as 
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The  resulting  equation  for  A  is  given  by 
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z'*' 

Here,  is  the  absolute  group  velocity  in  a  stationary 

reference  frame,  given  by 


-  C *  -f-  u.  (y.z.C) 

A#  Q  X"  /V  f  /v^ 

At  this  stage,  no  assumption  has  been  made  concerning  the 
orientation  of  k  or  U  with  respect  to  stationary  Cartesian 
axes  x  and  y.  The  terms  (I)  have  been  given  previously  in 
the  present  form  (but  for  one  space  dimension)  by 
Tur pin ( 198 1 ) ,  Benmoussa ( 198 3 )  and  Turpin,  Benmoussa  and 
Mei(1983),  who  have  studied  the  effect  of  following  and 
adverse  currents  on  the  shoaling  of  Stokes  wave  packets  (in 
the  form  of  solitons)  in  one  dimension.  The  steady-state, 
current-free  counterpart  of  the  term  (I)  is  directly  related 
to  the  spatially  varying  shoaling  term  in  the  parabolic 
approximation,  as  will  be  seen  below.  The  terms  (II), 
involving  derivatives  of  J  ,  represent  further  corrections 

of  O(o)  due  to  the  presence  of  a  varying  current.  The 

,  / 

presence  of  the  second  unknown  *pL  requires  a  second 

equation  for  closure;  this  is  given  by  the  forced  wave 

,  / 

equation  (3.3.18)  for  ^  . 
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A  number  of  virtually  independent  simplifying 
assumptions  may  be  applied  to  the  pair  of  coupled  equations 
(4.3.5)  and  (4.3.7),  and  we  illustrate  several  choices  here 
in  order  to  relate  the  results  of  the  present  study  to  those 
of  other  investigators.  The  first  thorough  study  of  Stokes 
waves  in  a  domain  with  ’’fast"  modulations  (  a 0(c)  )  was 
made  by  Chu  and  Mei (1970a ,b) ,  who  investigated  essentially 
the  same  problem  as  studied  here  but  neglected  the  presence 
of  a  large  current.  Dropping  the  terms  involving  the  0(1) 
mean  current  JJ  from  (4.3.5)  and  (4.3.7)  leads  to  the 
simplified  set 
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where  F1  is  a  simplified  form  of  F  obtained  by  setting  ^ 
and  £7=0.  Equa t ions ( 4 . 3. 8-9 )  contain  the  same  information 
as  the  results  of  Chu  and  Mei ,  who  left  their  equations  in 
the  form  of  a  set  of  conservation  equations. 


A  mild  slope  (  o  ■'d*' )  form  of  the  current-free  model 
may  be  constructed  by  neglecting  terms  of  0(  )  and  higher 
in  (4. 3.8-9)  after  the  shift  in  scaling.  We  remark  that  the 
modulations  of  A  must  be  left  at  0(fc),  since  nonlinear 
instabilities  such  as  the  Ben j amin-Fei r  instability  may  lead 
to  fast  amplitude  modulation  in  an  otherwise  slowly  varying 
domain  (see  Benjamin,  1967,  and  Benjamin  and  Fair,  1967). 
The  coupled  equations  then  further  reduce  to 
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with  (4.3.9)  altered  to 
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The  three  models  mentioned  so  far  may  in  principle 
be  solved  as  an  initial  boundary  value  problem  in  an 


1  -  iV.fMV  .1  1  mi  - 


70 


arbitrary  domain;  however,  the  determination  of  the 
direction  of  k  at  a  local  point  in  the  domain  is  problematic 
since  the  phase  surface  ^'(x,t)  is  not  being  determined. 
For  problems  in  domains  where  a  principle  direction  of 
propagation  may  be  defined,  it  is  convenient  at  this  stage 
of  the  derivation  to  specify  that  the  phase  surface  is  one 
dimensional  in  space.  Choosing  x  as  the  principle  direction 
of  propagation  and  defining 


=  kU‘)<Oi'  -  Cot 


leads  to  the  simplifications 
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The  amplitude  is  then  allowed  to  vary  in  the  y-direction  to 
account  for  modulations  and  small  deviations  from  the 
principle  propagation  direction.  The  ability  to  model  the 
deviation  in  travel  direction  will  be  crucial  to  the  results 
of  Chapter  6.  Inserting  the  simplifications  (4.3.12-13) 
into  (4.3.10-11),  and  noting  that 
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in  the  absence  of  a  current,  reduces  the  mild  slope 


formulation  to 
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which  is  the  extension  of  the  equation  of  Djordjevic  and 
Redekopp  (1978)  to  two  space  dimensions.  Note  that  the  term 
(CC^ A„ ),.  in  (4.3.15)  can  be  expanded  as 

=  C(Ta  +  ^CC^^Aty 


The  second  terra  is  formally  of  0(  o  )  in  the  mild  slope 
approximation,  and  could  be  dropped.  However,  we  will 
retain  it  for  the  sake  of  consistency  with  previously 
derived  models,  as  will  be  demonstrated  in  section  4.4. 
Finally,  for  water  of  constant  depth,  (4.3.15-16)  take  the 
following  form  after  substituting  for  time  derivatives  of  A 
using  (4.3.3): 
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Equations  (4.3.17-4.3.18)  are  equivalent  to  the  coupled 


equations 

d  e  r  i  v  ed 

by  Davey 

and 

Stewar  tson 

(1974 ) 

and 

Yue (1980) . 

It  is 

generally 

impo 

ssible  to 

el im inate 

,  t 

between  (4 

.  3. 17  )  and 

(4.3.18), 

and 
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so  must 

be  solved 

for 

in  addition  to  A.  The  exception  to  this  rule  is  in  the  case 
of  a  strictly  one  dimensional  envelope  A  oriented  colinearly 
or  at  some  angle  to  the  propagation  direction  x.  The 
colinear  case  corresponds  to  the  strictly  one-  dimensional 
problem  described  originally  by  Hasimoto  and  Ono(1972).  In 
the  case  of  the  envelope  function  A  making  an  oblique  angle 
with  the  direction  of  wave  propagation,  solutions  in  the 
form  of  oblique  envelope  solitons  and  groups  have  been  found 
by  Saffman  and  Yuen(1978  )  and  Hui  and  Ham il ton ( 1 97 9 )  ,  for 
the  case  of  deep  water,  and  by  Kirby  and  Da  1 r ymple ( 1 98 3c) 
for  the  case  of  intermediate  water  depth. 


For  cases  where  the  0(1)  current  is  retained,  the 
same  series  of  gradual  reductions  may  be  applied.  We 
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present  here,  in  order  to  show  consistency  with  previous 
investigators'  results,  the  equations  resulting  from  the 
assumption  of  a  mild  slope  and  propagation  in  one  spatial 
direction,  including  the  effect  of  a  following  or  opposing 
current . 


Dropping  higher  order  terms  (  II  )  and  the 
y-dependence  from  (4.3.5)  leads  to  the  one-dimension 
evolution  equation 


where  we  have  used  the  definition  (4.3.14)  for  7^  after 
replacing  «;  by  J"  ,  and  have  neglected  dissipation.  The 
second  derivative  terms  may  be  combined  using  the  lowest 
order  energy  equation.  Restricting  attention  to  a  slowly 
varying  (in  time)  current  field  then  reduces  (4.3.19)  to 


x-  (At  *  Cy  Ay  )  +  Ajr  fcpLj  A 

-  crklb!A]lA  - 

1  l  *  2rc; 


- 

17  C' 


In  f :  -  kL  AamI  k#  )  ?  A 


iL 


eoih 


'  Z’-L 


xx 


—  QA'lA  «£> 
ru  ~ot  J 
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The  corresponding  forced  wave  equation  for  p  is  given  by 

4't£  +  2a  4'xt  +  <4  4*  +  ^l4'x  )*  -  ^  ^  4X  )*  = 

2 

_ ok '4  I2-  <3  fttklM  \  ,  !a, 2) 

"  zzzijfi;  lAit  -rlT^rar/x  l#  1  ;x  6a  zO 

Equations  (4.3.20-21)  contain  the  same  information  as  the 
evolution  equation  of  Turpin,  Benmoussa  and  Mei  (1983),  who 
derived  a  single  equation  in  A  directly.  Defining  the 
moving  coordinate 


Tr  «  * 


f  I  -  ~t  1 
i  '  <=J- 


and  the  slow  coordinate 


Xi.  =  *'* 


(4.3.20-21)  may  be  combined  to  give 


A.  *  jr_  fCu\  A 


I—  (Si  )  f  i-  sk 


1  2rc,jc;jl:  %  C'd'ti  i 

+A.  <r£Zf D  - £-i—  p¥=3. —  +/)  — i — y—  l!-AlZA  * 

+  =0  6^-h) 

where  Q  is  a  real  function  involving  the  integration 
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constant  of  .  In  the  case  of  an  isolated  wave  packet 

l  f 

where  ^  0  for  x-v±-«  ,  Q  may  be  set  equal  to  zero 

Equation  (4.3.21)  is  equivalent  to  the  result  of  Turpin 
Benmoussa  and  Mei(1983). 
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4.4.  Parabolic  Equations 


Under  a  restrictive  set  of  assumptions,  the  elliptic 
model  of  section  4.2  and  the  hyperbolic  models  of  section 
4.3  may  each  be  reduced  to  parabolic,  initial  boundary  value 
problems  for  waves  propagating  in  or  close  to  a 
pre-speci f ied  direction.  The  first  assumption,  already 
applied  in  section  4.2,  is  that  the  wave  field  is  purely 
harmonic  in  time.  The  resulting  neglect  of  the  time 
dependence  of  the  amplitude  A  bars  the  consideration  of  wave 
fields  which  are  groupy  in  nature  as  well  as  precluding  the 
study  of  the  onset  of  instabilities.  The  parabolic 
equations  which  follow  are  thus  applicable  only  to  the  study 
of  plane  waves . 


The  second  major  assumption  involves  the 

specification  of  one-  directional  propagation.  This 

assumption  has  already  been  made  with  regard  to  the 

equations  of  section  4.3,  and  the  parabolic  equation  then 

follows  directly  from  these  results  after  neglecting  time 

dependence  and  introducing  a  shift  in  scaling.  The 

application  of  a  splitting  matrix  to  the  elliptic  model  in 

section  4.2  leads  to  a  set  of  coupled  equations  for  the 

\ 

forward  and  back-scattered  components  of  T'(  ,  as  in  the 

derivation  of  Radder  (1979).  In  this  section,  we  will  use 
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the  method  employed  by  Booij(1981)  to  obtain  a  single 

'T  + 

equation  in  the  forward  scattered  component  p,  which 

/y  “ 

neglects  the  coupling  with  the  back-scattered  component  r'(  . 
The  result  will  be  to  extent  the  parabolic  approximation  of 
Booij  to  include  the  effect  of  nonlinearity. 


As  a  first  example,  the  derivation  of  a  parabolic 
equation  from  the  mild-slope,  time-dependent  equation 
(4.3.15)  is  illustrated.  A  shift  in  scaling  is  introduced 

by  specifying  that  the  x-derivative  of  A  be  0(i5), 

« 

consistant  with  the  assumption  that  the  variation  of  the 


pnase 


accounts  for  all  the  fast  variation  in  the 


direction  of  propagation.  In  contrast,  the  y-derivative  of 

A  will  remain  at  0(:  ) ,  on  the  assumption  that  waves  turned 

at  a  small  angle  to  the  x-direction  will  exhibit  slow 

phase-like  variation  in  the  y  direction.  The  term  {CC^A^)^ 

will  be  kept  in  its  entirety;  this  choice  leads  to  the 

retention  of  the  formally  small  componant  (CC(,)„AM,  but  we 

note  that,  for  waves  turned  at  large  angles  to  the  principle 

propagation  direction,  A^  should  approach  0(e)  »  in  which 

case  (CC,j)^A^  is  comparable  in  order  to  the  remaining  terms. 

.  i  i 

We  also  note  that  T.  is  formally  of  0(c  )  in  this 
approximation.  These  scaling  assumptions  and  the  neglect  of 
time  dependence  reduce  (4.3.15)  to  the  parabolic  equation 


78 


+4(C,\A 


+  l  (CCd  Aca  )u  —  coJ?  D  l A  i  A  *  j-y.  A  =  o 

lu>  *  *  °  2.  x 


It  is  helpful  in  applications  of  parabolic  equations 
to  define  the  amplitude  with  reference  to  a  fixed  wavenumber 
kD  based  on  some  characteristic  of  the  physical  domain  of 
interest  (  usually  the  initial  conditions).  The 
introduction  of  k0  facilitates  the  reconstruction  of  the 
instantaneous  water  surface  since  all  phase  information 
not  contained  in  k0x  is  then  absorbed  by  the  amplitude 
function  A.  We  therefore  let 


x 

,  x  -  frdx) 

A  c  A  «f 


(v.y.O 


A'  is  then  the  amplitude  of  a  wave  given  by 

,  x.  i  x  —  ^  ) 

9,  -•  A  (Xj1)')  € 


(y.v.j) 


tx> 


Substituting  (4.4.2)  into  (4.4.1)  leads  to  the  revised 
equation 

A  A  x  +  L  £-  fct)C^  A  +  Cp  ^  A  -*•  _i-  A ^  -V 

-  0»kzVlfl\lA'  +  U&A'  =  o  6^*/} 

T  2- 


After  dropping  the  dissipation  term,  (4.4.1)  and  (4.4.4)  are 


^ nAA  in 


Util 


equivalent  to  equations  (2.18)  and  (2.22)  of  Kirby  and 
Dalrymple ( 198 3b) ,  which  were  derived  using  a  multiple  scale 
perturbation  technique.  Each  model  may  be  further  reduced 
to  an  equation  for  pure  diffraction  by  setting  k=k0=constant 
and  C,  C(j  =  constant.  Noting  that  A'=  A  in  this  case, 
(4.4.1)  and  (4.4.4)  both  reduce  to 

2+k0Ay  +  /4  -  uk)  01a11A  + wA  'O  (v.v.r) 

which  is  equivalent  to  the  diffraction  equation  of  Yue  and 
Mei(1980)  after  neglecting  the  dissipation  term. 

Kirby  and  Dalrymple  (1983b)  have  shown  that  the 
linearized  form  of  (4.4.1)  or  (4.4.4)  is  consistant  with  the 
parabolic  equation  of  Radder(1979)  for  the  forward  scattered 

I 

wave  ^  after  neglecting  dissipation.  This  consistancy  will 
be  demonstrated  again  subsequent  to  the  application  of  a 
splitting  method  to  the  elliptic  model  (4.2.3).  The 
equivalence  of  the  results  of  the  splitting  method  to  the 
t  ime- independent  form  of  the  nonlinear  Schrodinger  equation 
after  introducing  the  slower  scaling  in  x  allows  us  to  draw 
an  important  conclusion  before  proceeding  to  the  derivation 
of  a  general  model  including  currents.  If  we  neglect  the 
mean  current  U  and  the  time  dependency  in  the  wave-action 
equation  (4.3.3),  and  assume  that  ,  (4.3.3)  reduces 
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to 

c,  Ay  ♦  x  c,xA  -o  (v.y.4) 

which  is  equivalent  to  the  leading  two  terms  in  (4.4.1).  It 
is  apparent  that  a  requirement  for  a  correct  parabolic 
approximation  is  that  the  governing  equation  must  contain 
the  correct  form  of  the  relation  for  wave  action 
conservation.  To  further  illustrate  this  point,  we  drop  the 
higher  order  terms  (II),  introduce  the  slow  scaling  in  x, 
and  neglect  time  dependence  in  the  general  nonlinear 
Schrbdinger  equation  (4.3.5)  to  arrive  at  the  form 


It  is  apparent  that  the  bracketed  term  is  the  equation  for 
wave  action  conservation  including  energy  dissipation,  and 
contains  the  essential  features  of  (4.3.3).  We  may  require 
that  the  results  of  a  splitting  method  applied  to  the 


r 

■% 
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elliptic  model  (4.2.3)  should  reproduce  this  information  in 
as  close  to  the  exact  form  as  possible. 

Booij(1981)  presented  a  splitting  method  and 

y* 

resulting  parabolic  equation  for  the  velocity  potential  <£( 

based  on  the  mild-slope,  linearized  equation  (4.2.4).  A 

quick  inspection  of  Booij's  resulting  parabolic  equation 

shows  that  it  cannot  reproduce  the  terms  (UAV  +  VA„  )  and 

x  3 

[(U/~)  +  (V/j")  ]  which  are  important  components  of  the 

X  4 

wave  action  conservation  equation.  We  use  instead  a 
modified  form  of  Booij's  derivation  which  leads  to  a  more 
nearly  correct  form  of  the  wave  action  conservation 
equation . 


Booij  showed  that  the  second  order  elliptic  equation 


can  be  split  exactly  into  equations  governing  forward  and 

i  -t 

backward  moving  (with  respect  to  x)  disturbances  Ti>i  and 
which  are  given  by 


=  ji 

17  ^  •«* 


(V.v.%) 
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l  “* 

Since  the  split  is  exact,  the  equation  governing  is  then 

discarded  since  there  is  no  coupling  between  the  two 
potentials.  The  assumed  exactness  of  the  split  of  the  | 

i 

elliptic  model  and  the  subsequent  neglect  of  certain  I 

derivative  terms  removes  the  ODtion  of  calculating  the  I 

t 

reflected  wave  by  an  iterative  procedure,  as  employed  by  Liu  | 

and  Tsay(1983)  to  calculate  reflection  from  isolated  fc 

{ 

submerged  obstacles  in  the  absence  of  a  current.  The  ^ 

V  '  1  * 

remainder  of  the  problem  centers  on  determining  u  and  the  - 

i  ** 

relation  between  anc^  the  wave  potential  .  Booij  . 

chose  the  relation 

i 

I 

5?.  to*")  I’ 

r 

*-•*  f 

where  <f  is  an  as  yet  unknown  operator  on  .  Substituting  ; 

(4.4.10)  into  (4.4.8),  expanding  and  neglecting  derivatives 
of  y  alone  then  leads  to  the  result 


The  elliptic  model  (4.2.3)  must  now  be  put  into  the  form 
(4.4.11)  as  closely  as  possible.  Letting  p=CC<^  and 
expanding  some  terms  in  their  spatial  components  allows 
(4.2.3)  to  be  written  as 


i 


(pi, A  npl^l*  +  kp<p,  +  (tez-ez ++cofa-u)- <tzrdIa, 


zS 


**<r4  -^Ux 

-r  2*  ^  JA-Vh  k,  =  o 


W'H, 


.  / 

In  the  following,  we  will  neglect  the  interaction  with  Vj,^ 
and  b^_.  Booij  chose  to  drop  terms  involving  squares  of  the 
mean  current  components  based  on  the  assumption  that  they 
are  likely  to  be  arithmetically  smaller  than  terms  involving 
p  except  in  the  vicinity  of  stopping  points.  This  choice 
leads  directly  to  Booij's  inability  to  derive  the  correct 
wave  action  equation  and  also  would  lead  to  difficulties  in 
the  derivation  to  follow.  Booij  then  arranged  (4.4.12)  in 
the  form 


3>  +  p"  (o%~ Ho, u.)Z  JL  )  <A  «o 

~  XX  r  -  x 


/?> 


where  M  is  given  by 
r~ 


(Here,  Booij's  derivation  is  extended  to 


include  the 


S4 


nonlinear  term).  (4.4.13)  is  then  compared  to  (4.4.11)  to 

complete  the  splitting  process.  We  note  that  the  inclusion 

of  the  term  2^*0^  in  the  splitting  procedure  leads  to  the 

neglect  of  the  term  LJA x  in  the  final  parabolic  equation, 

since  UA^  is  already  contained  in  the  term  .  We 

thus  follow  an  alternate  approach.  Assuming  that<  waves  are 

:  f  b  4_  x 

oriented  in  the  +x  direction  (phase  given  by  e  ’’  )  ,  the 

dispersion  relation  may  be  written  as 

cr  -  CO-  kU 

leading  to  the  result 

coz-crl-=  2 cjkli  -  [km)1' 

Equation  (4.4.12)  may  then  be  written  as 

Jx  +  r\p-'S)  £  +  =  o  >.^./7) 

/y 

where  M  T|  is  given  by 

C  1 

M  S,  =  [looktJL  +  -  crzk  D!aIZ  +sC7k/  j  J>, 

(v V.  tt) 
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In  the  following,  we  will  neglect  the  interaction  with  ^4^ 
and  b^.  Booij  chose  to  drop  terms  involving  squares  of  the 
mean  current  components  based  on  the  assumption  that  they 
are  likely  to  be  arithmetically  smaller  than  terms  involving 
p  except  in  the  vicinity  of  stopping  points.  This  choice 
leads  directly  to  Booij' s  inability  to  derive  the  correct 
wave  action  equation  and  also  would  lead  to  difficulties  in 
the  derivation  to  follow.  Booij  then  arranged  (4.4.12)  in 
the  form 

4>  xx  +  p"  (pK-Zu.u)^  +kx(if  I)^o  (M.12) 

fS* 

where  M  is  given  by 

^  \  /v 

M^>  •  sx  +/iuf7h )- I  +*rv  )  <p,  *■ 


(Here,  Booij 's  derivation  is  extended  to 


include  the 
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nonlinear  term).  (4.4.13)  is  then  compared  to  (4.4.11)  to 
complete  the  splitting  process.  We  note  that  the  inclusion 
of  the  term  in  the  splitting  procedure  leads  to  the 

neglect  of  the  term  UAX  in  the  final  parabolic  equation, 

rf 

since  UA^  is  already  contained  in  the  term  x  .  We 

thus  follow  an  alternate  approach.  Assuming  that^  waves  are 

.•  ffc 

oriented  in  the  +x  direction  (phase  given  by  e  '*■  )  ,  the 

dispersion  relation  may  be  written  as 


to  -  hU 


fa. its) 


leading  to  the  result 


u>z~crz*  2cj  kU.  -  fal/J) 


fa.i.ii) 


Equation  (4.4.12)  may  then  be  written  as 


fa.H.n) 


where  M is  given  by 


*  j2  laklL  +ACo(j?h  '&)  -  CTzkl  DM  I  Z  +/trh/]  <j>, 

-(wLl  )  +£ (p-vl)Ll  *&*>.<*■%$ 


fa.i./t) 


We  have  thus  imbedded  the  term  U  •£.  ,,  which  we  don't  want  to 

C* 

reduce  in  order,  in  M  £>  .  (4.4.17)  may  then  be  written  in 
the  form  of  (4.4.11); 

?.xx  +  ?,x  +  + 

Comparison  with  (4.4.11)  leads  immediately  to  the  results 


*  -  k  (\  + 

(y.i.io) 

?2/y  -  (p-  Ur) 

&Ut) 

y-  #?(i+ 

rtp-u*) 

(y.y.  22) 

? *  Cp-u)'(i+  k'(f-u.') )' 
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Booij  points  out  that  Radder 1 s (1979 )  approximation  is 
equivalent  to  the  immediate  choice 


while  Booij's  approximation  retains  the  full  expression  and 
then  approximates  the  resulting  parabolic  equation 


The  results  of  both  expansions  may  be  summarized  in  the 
general  equation 


? 


N 


where 


csiv  ar.-i  l-.  v  v 


and  where  Pi  -P(  ■ 1/2  for  both  approximations.  The 
approximation  with  P(  *0  then  should  correspond  to  the 
information  contained  in  the  nonlinear  Schr'odinger  equation. 
Equation  (4.4.29)  may  be  expanded  to  yield  the  expression 


+  x{kcP-ui)]j;  » 

*  A  k^Cp-U,1)  $+  +  ^  + 


&Y.  v) 


Then,  an  equation  for  the  complex  amplitude  is  obtained  by 


'%  T 

expanding  M and  making  the  substitution 

V  *  -*(&) 2  ^ 

which  results  in  the  equation 


(V.y.jO 


(c**  u)a*  *  (f)  v\  * 

-  +  H  A  +a<t£.z  DjAlzA  =  o  (y. y. ii) 


where  terms  containing  products  of  components  of  the  ambient 
current  have  been  dropped.  This  equation  may  then  be 


compared  with  (4.4.7),  and  we  note  that  errors  have  been 
introduced  through  the  factor  {*0/0')  multiplying  VA^  and 
the  appearance  of  U)  rather  than  <7*  in  the  y  derivative  term 
for  V  and  J" .  A  complicated  procedure  involving  a  double 
splitting  of  the  equation  in  x  and  y  shows  that  the  errors 
appearing  in  (4.4.33)  result  from  the  choice  of  Jk»k^x  and 
the  effect  of  this  choice  on  the  assumed  dispersion  relation 
(4.4.15).  However,  the  resulting  approximation  is  more 
complete  than  the  approximation  obtained  by  Booij  due  to  the 
inclusion  of  the  term  UA^  and  the  correct  form  of  the 
derivative  of  (C^  +U)/d".  Errors  resulting  from  the  use  of 
(4.4.33)  rather  than  the  more  strictly  correct  form  (4.4.7) 
will  be  investigated  in  section  6.9.  The  approximate 
equation  (4.4,34)  is  further  modified  by  the  shift  to  a 
reference  phase  function  e ,  resulting  in  the  equation 

+  -&{<** A'*\*  +  + 


+  /L  DlAlxA‘ '  *=  O 

2. 


(V.y.3  % 


where  A'  is  related  to  by  (4.4.2).  In  the  absence  of  a 
current  (£,  (4.4.34)  reduces  identically  to  (4.4.4). 


The  equations  (4.4.4)  and  (4.4.34)  form  the  basis 
for  the  computational  schemes  used  in  most  of  the  examples 
presented  in  Chapter  6.  In  several  examples  where  currents 
are  neglected,  use  is  also  made  of  the  more  accurate 
approximation  with  P(  *1/4 ,  P^*3/4.  The  derivation  of  the 
equation  governing  the  complex  amplitude  A'  for  this  case 
follows  along  similar  lines;  however  several  simplifications 
are  introduced.  In  particular,  the  approximation 


2.  f 

>x  l 


mlVj 

CC+  J 


(v.r.Jf) 


is  made,  consistent  with  the  neglect  of  amplitude 
modulations  in  the  nonlinear  term.  Further,  we  introduce 


Jl  £  J.O-W  |  ^  -  o-kv^,  * 


consistant  with  the  assumption  of  a  small  coefficient  for 
the  damping  term.  Substituting  (4.4.32),  (4.4.35)  and 

(4.4.36)  into  (4.4.30)  with  Pj  *1/4  and  P^  *3/4  and 
neglecting  currents  leads  to  the  governing  equation 

/C  Ay  (t-h)C^A  ^  A  +  A ^  + 

* ± 5. Ccc^A,.  )  x  a  )  - 

tstk  ^  ^  uki  *  2klcc*  1 
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where  we  have  also  neglected  some  small  derivative  terms 
multiplying  the  nonlinear  and  damping  terms.  Equation 
(4.4.4)  is  recovered  from  (4.4.37)  by  setting  P,  «0  and  Pt 
■1/2,  consistant  with  Radder's  splitting  scheme.  The 
numerical  schemes  used  to  solve  (4.4.4),  (4.4.34)  and 

(4.4.37)  are  discussed  in  section  6.2. 


CHAPTER  5.  COMPARISON  OF  THE  MILO  AND  MODERATE 


SLOPE  APPROXIMATIONS  FOR  LINEAR  WAVES 

5.1.  Introduction 

For  most  applications  to  water  wave  propagation 
problems,  the  mild-slope  form  of  .the  models •  derived  in 
Chapter  3  should  be  adequate  for  a  reasonably  accurate 
representation  of  the  wave  field.  In  particular,  ocean 
bottoms  consisting  of  loose  material  typically  have  bottom 
slopes  on  the  order  of  0.001-  0.01,  and  retention  of  terms 
in  (Vfch)1  and  ^h  is  clearly  unnecessary.  Tidal  currents 
interacting  with  these  smooth  topographies  would  also  be 
expected  to  be  very  slowly  varying  in  the  horizontal 
directions,  thus  alleviating  the  need  to  retain  terms  to 
0(S  )  covering  variations  of  CT  and  k.  The  mild-slope 
approximation,  given  in  linear  form  by  Booij  (1981)  and  by 
(3.4.1)  and  in  nonlinear  form  by  (3.5.9),  should  then  be 
adequate  for  the  modelling  of  waves  in  most  physically 

realistic  domains.  The  principle  exception  to  this 

'y  rr 

conclusion  is  that  the  term  involving  in  F  ^  is 

required  for  the  study  of  unsteady  nonlinear  waves  in  order 


92 


to  avoid  tha  types  of  singularities  in  amplitude  observed  by 
Lighthill  (1965, 19€7)  in  his  study  of  Whitham's  conservation 
laws;  see  also  Chu  and  Mai  (1970a, 1971) .  This  requirement 
is  consistant  with  the  fact  that  the  amplitude  may  become  a 
fast  varying  quantity  even  in  a  very  slowly  varying  or 
homogeneous  domain  as  wave  instabilities  develop  and  wave 
groups  form.  The  connection  of  the  term  in  p, to  the 
higher  order  dispersive  terms  in  the  NLS  was  discussed  in 
Chapter  4. 


Exceptions  to  the  assumptions  of  a  globally 
slowly-varying  domain  can  arrise  in  several  situations. 
Relatively  large  bottom  slopes  may  exist  over  short  physical 
expanses  on  rocky  bottoms  or  in  the  vicinity  of  naturally  or 
artificially  maintained  channels,  and  tidal  currents 
interacting  with  these  relatively  rapidly  varying  bottoms 
will  themselves  exhibit  the  same  rapid  variations  provided 
that  the  changes  in  depth  occupy  a  significant  proportion  of 
the  total  depth.  In  addition,  ambient  currents  in  an 
otherwise  slowly  varying  domain  may  exhibit  fast  variation 
due  to  specialized  circumstances,  such  as  the  flow  fields 
characteristic  of  tidal  ebb  currents  from  inlets  and  rivers, 
rip  currents,  and  discharge  plumes  from  outfall  structures. 
Such  concentrated  flows,  which  often  assume  the  form  of 
turbulent  entraining  jets,  violate  the  assumptions  of  the 
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present  theory  in  two  important  ways.  First,  the  small  and 
intermediate  scale  turbulence  of  the  ambient  flow  field 
violates  the  assumption  of  irrotational ity  on  a  local  scale 
of  the  order  of  a  wavelength  or  smaller.  Savitsky (1970) 
studied  the  refraction  of  waves  in  the  presence  of  a 
turbulent  flow  with  mean  current  gradient,  and  concluded 
that  the  effect  of  the  turbulence  was  small  in  comparison  to 
the  interaction  with  the  mean  characteristics  of  the  flow. 
This  result  is  encouraging;  however,  the  mean  flows 
themselves  are  generally  rotational  over  a  large  scale.  It 
is  known  (see,  for  example,  Mei(1982),  section  (3.6)  )  that, 
for  current  variations  of  0(  €*■  )  ,  the  refraction 
approximation  for  waves  in  the  presence  of  rotational 
currents  is  unaltered  from  the  strictly  irrotational  case. 
No  such  results  are  currently  available  for  the  faster 
variations  implied  in  this  study. 


In  this  chapter,  we  study  several  cases  involving 
relatively  steep  bottoms  with  the  intent  of  determining  the 
range  of  significance,  if  any,  of  the  0(£L)  terms  added  to 
the  mild  slope  form  in  Chapter  3. 
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5.2.  Linear  Edge  Waves  / 

As  a  first  example,  we  investigate  the  effect  of  the 
steep  bottom  terms  on  numerical  predictions  of  the  longshore 
wavelength  of  edge  waves.  We  consider  the  topography 

In  (x)  =  'Wt  X  (XZ.l) 

where  m  is  the  beach  slope  and  x  is  directed  offshore. 
Edgewaves  are  trapped  modes  which  propagate  alongshore  in 
the  +  or  -  y  direction  (  or  both,  in  the  case  of  standing 
waves  )  and  are  trapped  between  the  shoreline  and  an 
offshore  turning  point  by  refraction.  For  water  of 
intermediate  depth,  CJrsell  (1952)  has  provided  an  exact 
solution  for  the  discrete  spectrum  of  the  linearized 
solution;  the  dispersion  relation  is  given  by 

CJ1  ,  sin  [  (  2 -h  +  /  )  +An  * >•  ]  (S.  2.2.) 

in  the  absence  of  currents,  where  ?.  is  the  longshore 
wavenumber  of  the  edge  wave  and  n  is  the  mode  number.  A 
family  of  normalized  surface  profiles  for  n=0-5  is  shown  in 
Figure  5.1.  The  trapping  mechanism  requires  that  the 
longshore  wavelength  be  shorter  than  the  deepwater 
wavelength  for  a  wave  of  the  same  frequency,  or 


V. 
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3  ^  Jal  =  fr.  Cs.  2.j) 

“V 

where  k0  is  the  deepwater  wavenumber.  For  each  given  mode 
number  n,  (5.2.3)  then  requires  that 

SlU  £  (I'K-h  0  f Ik  ]  <  1  (S’.Z.i) 

or 

'H  ■<  ^AA/‘f  7*  ?  (T.  2.r) 

4  1727^7)  J 

Condition  (5.2.5)  must  be  met  for  a  mode  n  wave  to  exist  on 
a  given  beach.  Interestingly,  (5.2.5)  is  satisfied  for  all 
finite  values  of  m  when  n*0. 


% 


In  the  present  application,  we  study  the  prediction  ip 

of  /'.  by  reducing  the  linear  wave  equation  to  an  ordinary 
differential  equation  on  the  domain  7\  is  then 

related  to  the  eigenvalues  of  the  reduced  system.  After 

, 

„  n 

neglecting  currents  and  dissipation,  and  assuming  purely 
harmonic  motion,  (3.3.8)  reduces  to 

£2 

]  +  klcC}$  -  ^C<=L+  F-)?,  -o  (S.II-)  :i 
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The  coefficient  F  involves  the  expression  k.^(  )  in  each 
term.  If  we  interpret  k  to  be  oriented  in  the  direction  of 
A  alongshore,  it  is  apparent  that  each  term  in  F  is  zero, 
and  (5.2.6)  reduces  further  to  the  form 


The  mild  slope  equation  is  recovered  by  dropping  oi  .  Smith 
and  Spr inks (1975)  used  the  mild  slope  equation  to  calculate 
A  on  a  plane  slope,  and  Kirby,  Dalrymple  and  Liu  (1981) 
have  used  the  mild  slope  equation  in  a  finite  difference 
form  to  study  edge  waves  on  arbitrary  bottom  profiles.  In 
the  present  case,  we  will  use  (5.2.7)  to  determine  the 
effect  of  the  term  oi  on  the  prediction  of  values  of  A  , 
using  a  numerical  method  similar  to  that  of  Kirby,  Dalrymple 
and  Liu. 


r  — 

First,  the  wave  ^  is  assumed  to  be  of  the  form 
ry  . 

1  /  \  1  /_  \  »  o  i 


4s£m>  =  £,60c 


(r.2.r) 


Substituting  (5.2.8)  into  (5.2.7)  yields  an  eigenvalue 
problem 


ecu: 
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Boundary  conditions  on 


* 


are  given  by 


r 

b  ouulsJL  as  X  ■*»  0 

' 

L  0  AS  X  ■*  oo 


(S.  2.  lo) 


The  shoreline  boundary  condition  is  computationally  awkward, 
and  we  introduce  a  new  dependent  variable 

?  =  v?, 


such  that  (5.2.9}  and  (5.2.10)  become 

<***l?»  +  { -  2CSX 1  ?*  + 

+  -tCccJx  *  xYfelcC}-^rf)?f  *  **«*? 

Cr.  i.  a ) 


f  =  o  j  (r.z.it) 

Nondimensional  variables  are  introduced  according  to 

Oi.x,  '/I,' /*)-A(ir.V/V/S) 

401 


)  *  -S.  f  ^ )  c<5 ) 


d  3  d. 


where  ”  denotes  dimensionless  quantities.  Dropping  **'s  and 
defining  p»CC_  reduces  (5.2.11)  to 

pxlfxl<  *  CxlPx  -^xpl  ?x  +  i1?-  XP* 

=  llxlp  £  (s. a; 

The  dependent  variable  x  is  then  stretched  according  to 

X  «  e  -  1  (f.2.li 

in  order  to  provide  greater  detail  near  the  shoreline, 
yielding  the  final  equation 

POn-l)>„  - >[(«*•-!)%  2eW. 

3p]?, 

+  y  lczii(ex'-i)1  (*-*)]$  = 
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Equation  (5.2.15)  is  solved  on  a  truncated  domain 
0£z£zmax,  where  zmax  is  cnosen  such  that  the  wave  profile 
for  the  highest  desired  edge  wave  mode  may  decay  in  an 
unconstrained  fasmon,  with  ^(zmax) *0.  Tne  numerical  scheme 
employs  central  differences  on  eveniy  spaced  points  zx. 

Results  for  the  dimensionless  inverse  wavenumber 
U)/q as  a  function  of  beach  slope  m  are  plotted  for  modes 
0/1,  and  2  in  Figure  5.2.  Results  for  both  the  mild-slope 
(neglecting  oi. )  and  the  moderate  slope  cases  are  virtually 
indistinguishable  from  the  exact  solution  for  modes  1  and  2. 
For  the  mode  0  wave,  results  were  calculated  up  to  a  beach 
slope  m*2.Q.  sotn  the  moderate  and  mild  slope  results  are 
seen  to  diverge  from  the  exact  solution  witn  increasing 
slope,  witn  the  mild-slope  model  overpredicting  the 
wavenumoer  and  tne  moderate  slope  model  underpredicting  it. 
The  moderate  slope  results  are  seen  to  diverge  from  the 
exact  solution  earlier  and  faster  chan  the  mild  slope 
results.  Table  5.1  presents  a  comparison  of  calculated 
values  of  Lc  /'g  /*  from  each  model  witn  tne  exact  results. 

A  possible  explanation  for  the  faster  divergence  of 
the  moderate  slope  results  from  the  exact  solution  in  the 


range  m>.5  is  that  tne  large  slope  values  invalidate  the 
expansion  procedure  used  to  ootain  the  term  CL,  in  whicn  the 
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slope 

m 


.05 
.1 
.2 
.5 
0.75 


1.25 
1.50 
1.75 
2. 


w2/gi 

Ursell 

(1952) 

w2/gl 

mild 

slope 

% 

Error 

mild 

slope 

w2/gX 

moderate 

slope 

% 

Error 

moderate 

slope 

.01000 

.01023 

2.31 

.01023 

2.31 

.04994 

.05078 

1.69 

.05080 

1.72 

.09950 

.10070 

1.21 

.10084 

1.34 

.19612 

.19739 

0.65 

.19849 

1.21 

.4472 

.44583 

-0.31 

.45977 

2.81 

.6000 

.59217 

-1.30 

.62668 

4.45 

.7071 

.69061 

-2.33 

.74499 

5.36 

.7809 

.75717 

-3.04 

.82196 

5.26 

.8321 

.80418 

-3.36 

.88036 

5.80 

.8682 

.83872 

-3.40 

.91760 

5.69 

.8944 

.86490 

-3.30 

.94400 

5.52 

slope  is  assumed  to  be  a  small  parameter.  Both  models 
exhibit  moderately  large  errors  for  small  slope  values;  this 
result  is  probably  due  to  numerical  discretization  errors. 
The  error  relative  to  the  exact  solution  of  Ursell(1952)  is 
shown  in  Figure  5.3  for  the  choices  of  N=50  and  90  and 
zH**a I*  further  increase  in  N  reduced  the  error  for  small 
values  of  bottom  slope  only  slightly,  indicating  slow 
convergence  of  the  numerical  solution. 

The  numerically  calculated  values  of  X  differ 
slightly  in  the  range  0.03<m<0.2.  For  higher  values  of  m, 
the  two  solutions  diverge  rapidly.  Due  to  the  presence  of 
the  numerical  error  at  small  values  of  m,  it  was  impossible 
to  determine  which  model  (mild  or  moderate  slope)  exhibits 
more  accurate  behavior  for  intermediate  values  of  the  bottom 


5.3.  Reflection  from  Underwater  Topography 

In  this  section,  we  study  the  reflection  and 
transmission  characteristics  for  one  case  of  plane  waves 
incident  on  one-  dimensional  bottom  variations.  The 
topography  will  be  taken  to  vary  in  the  x-direction,  and  is 
assumed  to  be  uniform  in  the  y-direction.  Plane  waves  are 
incident  from  x»-  to  at  an  angle  $  to  the  x  axis.  The 
variations  of  the  topography  are  considered  to  be  of  bounded 
extent  horizontally;  however,  the  depths  at  ±00  are  allowed 
to  differ.  The  physical  domain  is  described  schematically 
in  Figure  5.4.  In  the  limit  of  topographic  variations 
characterized  by  vertical  walls  separating  regions  of 
constant  depth,  the  general  case  of  differing  deptns  at  ±00 
has  oeen  studied  by  Lassiter  (1972)  for  normal  wave 
incidence  (d»0)  and  by  Kirby  and  Dalrympie  (1983a)  for 
arbitrary  angles  of  incidence,  using  metnods  based  on  the 
matching  of  eigenfunction  expansions  along  the  vertical 
boundaries  of  subdomains  of  the  fluid  region.  Kirby  and 
Dalrympie  also  included  results  obtained  using  the  boundary 
integral  method  (BIM)  of  Raichlen  and  Lee (1978)  as  a  test  of 
the  principle  results,  and  extended  the  method  to  include 
oblique  wave  incidence. 

A  wave  propagation  method  based  on  the  mild  slope 


equation  (3.4.1)  has  been  applied  by  Booij(1981)  to  study 
the  reflection  and  transmission  characteristics  for  waves  in 
regions  characterised  by  gentle  slopes.  In  particular, 
B0013  studied  the  reflection  of  waves  normally  incident  on 
an  underwater  slope  joining  two  regions  of  constant  depth 
(with  waves  propagating  from  deep  to  shallow  water  and  being 
partially  reflected) ,  and  also  studied  oblique  incidence  on 
a  submerged  trench.  Of  particular  interest  in  this 
connection  is  the  case  where  a  caustic  forms  on  the  upwave 
trench  boundary.  If  the  trench  width  is  narrow  enougn  to 
allow  the  exponentially  decaying  wave  amplitude  in  the 
geometric  shadow  to  reach  the  second  caustic  line  on  the 
downwave  boundary,  then  significant  wave  energy  may  be 
transmitted  past  the  finite-width  shadow  boundary. 

We  wish  to  determine  the  importance  (or  the  ability 
to  improve  results  )  of  the  higher  order  terms  in  S  added 
to  the  mild  slope  formulation  in  Chapter  3.  The  present 
formulation  allows  for  unambiguous  specification  of  the  term 
Oi,  which  depends  only  on  the  bottom  slope  and  wave 
frequency.  However,  the  0(5*)  term  F  depends  on  the 
specification  of  a  wavenumber  vector  it  and  the  gradient  of 
the  amplitude  A,  and  these  quantities  are  not  well  defined 
in  a  partial  standing  wave  system.  The  assumptions  used  in 
evaluating  F  are  discussed  below.  We  neglect  dissipation 
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and  restrict  our  attention  to  linear  plane  waves  with  steady 
amplitude  in  the  absence  of  currents.  For  tnis  case,  the 
linear  model  (3.3.8)  can  be  written  as  (after  assuming 
harmonic  motion) 


vh-  +  s°  Cr-3-0 


The  restriction  to  one-dimensional  topography  allows 
for  some  major  simplifications  to  the  general  governing 
equation.  The  definition  of  the  wavenumber  vector  jc  leads 
to  the  result  that 

-O  iS.i.l) 

Then,  writing  k  in  its  components  fc*{kcosd  ,  ksin&}  leads 
to  the  result 

(  H  Si»  *  0 

or 


kM  s  w9(x)a  nyt  *  covshvi  (r.3.l) 


wnich  is  a  statement  of  Snell's  law 


Thus,  for  a  given 
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distribution  of  depth,  the  local  wave  angle  9  (x)  is 
uniquely  -determined.  The  wavenumber  component  X  in  the 
x-direction  can  then  be  determined  according  to 

JU),  (kh-^)‘A  frit) 

For  waves  propagating  from  shallow  to  deep  water,  it  is 
possible  for  m  to  be  >  k(x);  according  to  (5.3.4),  X  then 
takes  on  imaginary  values,  and  no  physical  significance  can 
be  attached  to  the  local  value  of  cos  9  (x)  .  The  point  =0 
is  a  turning  point  of  the  differential  equation  and 
represents  a  caustic  line  extending  in  the  y-direction 
through  the  point.  Waves  in  the  region  of  JL  imaginary  are 
constrained  to  travel  parallel  to  the  caustic,  with  the  wave 
profile  decaying  exponentially  away  from  the  turning  point. 

From  our  one  dimensional  assumption, 

Cr.  ir) 

and  the  elliptic  model  is  reduced  to  the  O.D.E.  : 

{ cc \  £  x  )x  +  X  CC}<£,  -  3  U  +  fH,  ■  O  (.s-.i.i. 

The  term  at  may  be  written  as 


T  ^ 
4*  *  i**  <p, 


vstGSMZXXMxt. 


irt.t'vneu 


S"  V  S  K-  *- 


n 


n 


0(  =  J,  ^  ^ l)x k„  +  S3k*  -  *  S,k,\  (r.3.7) 


The  term  F  is  evaluated  by  assuming  only  the 
presence  of  the  incident  wave  propagating  in  the  +x 
direction.  This  assumption  will  lead  to  errors  which 
increase  in  magnitude  with  increasing  reflection.  In  the 
absence  of  a  current,  F  may  be  written  as 


Fs  V1  -  V/4+/0 

'  ^1  +  ®/3  ■+  °<3  -2  3*j  ^ 


terms  ex',  ,  and  ^  may 


be  written  as 


Xhjh  5  =  Akjik'  •  fi.AA^/k'A 

ft  3  **  /■**'  ; 

?>' jtUKX/**  j 


A  V* 


*  w**  «'*  »"'*  j*  *  ■  *  •  ■** » *  ■  *  ' .*•  *, *  *„•'  \»  *  •  ■  »  *  •  *„«  "  •  *,»  *,»  *„•  * « *,*  *  •  *  *  “  •  *  ■  *  -  *  ■  ■  • 


/Sj  =  X  ( £A „  X  /k*A  .  /s  »  £(x*  X  A*v  (x.sa) 

The  gradient  of  A  may  be  further  evaluated  using  the  energy 
conservation  equation 


vh  -  (*£))  *o 

yielding 

(c.  XA*/k  )J(  »  O  (S.3.io) 

or 

Ay  «  -  a 

zc^x/a 


and  thus  A  may  be  eliminated  from  the  evaluation  of  F. 


Solutions  of  (5.3.6)  are  obtained  by  rewriting 
(5.3.6)  in  the  form 


4,  +  Ccc.X  Z  (£l~  « (x±£l  )4,  (ei/z) 

Tc^  T*  cc> 


and  then  writing  (5.3.12)  in  central  difference  form. 
Boundary  conditions  are  specified  at  x  positions  outside  the 


region  of  varying  depth.  In  this  case,  we  may  specify  the 
wavefield  on  the  incident  wave  (-x)  side  as 


r~> 

<j>  *  4>z  -  Cs.xa) 

where  is  the  incident  wave  of  unit  amplitude  and  <i>r  is 

a  reflected  wave  with  amplitude  K_ .  K_  is  then  the 

r  r 

reflection  coefficient.  The  transmitted  wave  at  x  downwave 
of  the  bottom  variation  is  given  by 

(s.in) 

Application  of  radiation  conditions  at  the  x  boundaries  then 
leads  to 


A.  A  (z?x 

A- A 


X  -=»  oo 

v  -ob  £<r.J./r) 


where  JLf  and  A  refer  to  values  of  at  h=h;(x*-o»)  and 
h»hx(x*+«e )  respectively.  Equations  (5.3.12)  and  (5.3.15) 
yield  a  tridiagonal  system  which  is  solved  by  a  standard 
double  elimination  method. 
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As  an  example  o£  reflection  from  an  isolated  deptn 
variation,  we  choose  a  bottom  consisting  of  a  smooth 
transition  from  a  aeptn  h( (x£0)  to  a  depth  h^(x£L),  with  the 
bottom  profile  in  the  region  0<x<L  given  by 

h(x)=  'lilh  OrM 

Z  Z 

where  s» (h( -h^) /L.  For  both  cases  studied,  we  take 
h^/hj-0.2,  with  waves  incident  from  the  depth  hj  .  The 
maximum  bottom  slope  occurs  at  x*L/2  and  has  the  value 
Tfs/2.  Results  were  calculated  for  values  of  s*l/3  and  2/3 
ana  for  angles  of  incidence  of  0*0*  and  45*.  Reflection 
coefficients  for  each  of  the  four  cases  are  calculated  by 
four  methods;  mild  slope  equation  (5.3.12)  neglecting  O f  and 
F,  moderate  slope  equation  (5.3.12)  retaining  only  the 
unamoiguous  term  Of ,  moderate  slope  equation  (5.3.12)  using 
of  and  F,  and  the  boundary  integral  method  (B1M)  described 
in  Appendix  F.  The  boundary  integral  method  provides  the 
theoretically  correct  answer  and  may  be  used  as  the  basis 
for  comparison  between  the  wave  propagation  methods  based  on 
(5.3.12) ,  although  the  metnod  described  in  Appendix  F 
neglects  local  curvature  of  the  bottom  boundary  and 
therefore  may  introduce  some  error.  The  discretised 


boundary  for  the  BIM  for  the  case  of  s*2/3  is  shown  in 
Figure  5.5. 


Results  for  the  reflection  coefficient  Kp  as  a 
function  of  kh }  for  the  incident  wave  are  shown  in  Figures 
5. 6-5. 9  for  the  four  cases  mentioned  above.  For  the  case  of 
s*l/3  (Figures  5.6  and  5.7),  the  mild  slope  equation 
provides  the  best  agreement  with  results  of  the  Bill  for  the 
45°  angle  of  incidence,  while,  for  normal  incidence,  the 
full  formulation  with  ei  and  F  agrees  more  closely  except  in 
the  range  of  kb1 <0.75,  where  Kp  is  large  and  errors 
associated  with  the  specification  of  F  become  important.  In 
tnis  range,  the  moderate  slope  equation  retaining  only  oL 
agrees  more  closely  with  the  BIM  results. 

For  the  case  of  s*2/3,  the  mild  slope  equation 
provides  the  best  estimate  of  the  reflection  coefficient  for 
normal  wave  incidence,  while  for  £*45*  the  moderate  slope 
equation  with  F  neglected  provides  the  best  agreement  over 
the  range  of  kh  values  tested.  Agreement  between  BIM  and 
the  moderate  slope  equation  with  ci  and  F  is  poor  for  both 
angles  of  incidence  in  tnis  case. 

The  results  of  the  cases  studied  here  indicate  that 
none  of  tne  approximate  methoas  are  accurate  in  terms  of 
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predicting  the  magnitude  or  variations  of  the  reflection 
coefficient  with  changing  incident  wavelength,  although  all 
of  the  methods  produce  qualitatively  similar  results 
approaching  the  limit  of  small  kh,  where  accurate  results 
may  be  obtained  simply  by  matching  surface  elevation  and 
mass  flux  across  a  discontinuity  in  depth  between  h=hj  and 
h^.  Of  the  models  described  above,  it  would  be  expected 
that  the  moderate  slope  model  with  F  retained  would  yield 
the  more  accurate  results  in  comparison  to  the  BIM,  although 
a  rational  scheme  for  specifying  F,  based  either  on  a 
formulation  for  standing  waves  or  an  iterative  scheme 
linking  the  incident  and  reflected  waves,  remains  to  be 
developed . 


CHAPTER  6.  SOME  EXAMPLES  OF  THE  COMBINED  REFRACTION- 


DIFFRACTION  OF  LINEAR  AND  STOKES  WAVES 


6.1.  Introduction 

In  this  Chapter,  we  use  the  time- independent 
parabolic  models  of  section  4.4  co  study  several  examples  of 
tne  combined  refraction  and  diffraction  of  linear  and  Stoxes 
waves  by  variable  dept n  and  currents.  Examples  include  a 
study  of  wave  focussing  by  a  submerged  snoal  and  diffraction 
of  waves  by  a  snore-attached  breakwater,  for  wnich  data  are 
available  for  comparison. 

Yue  and  Mei(1980)  have  snown  that  the  reflection  of 
a  wave  at  glancing  angle  of  incidence  on  a  vertical  wail 
produces  an  effect  known  as  the  "Mach"  stem,  where  the 
linear  result  of  an  incident  plane  wave  and  scattered 
reflected  wave,  leading  to  a  short-crested  pattern  in  tne 
far  field,  is  replaced  by  a  region  of  plane  waves  of  uniform 
amplitude  propagating  along  tne  wall.  Peregrine(1963) 
analysed  Yue  and  Mei's  parabolic  approximation  for 
diffraction  of  Stokes  waves  and  found  that  tne  equation  is 


121 


analogous  to  tne  nonlinear  equation  governing  tne 
development  of  an  undular  oore.  The  incident  wave  and  the 
wave  in  tne  stem  region  can  thus  be  interpreted  as  being 
con3ugate  wave  states  divided  by  a  jump.  Peregrine 
hypothesised  tnat  similar  effects  couia  be  expected  in 
situations  where  refraction  leads  to  ray-crossing  and 
overlapping  of  segments  of  tne  incident  wave.  Kirby  and 
Dalrymple { 1983b)  have  demonstrated  this  effect  in  the  case 
of  waves  focussed  by  a  submerged  snoal.  The  development  of 
a  conjugate  state  also  nas  implications  for  the  evolution  of 
wavefields  in  the  vicinity  of  straight  or  curved  caustics, 
wnere  the  development  of  a  jump  may  preclude  tne  evolution 
of  steady  state  solutions  for  the  wave  field  in  the  vicinity 
of  tne  caustic.  These  steady  state  solutions  can  be  found  in 
terms  of  the  second  Painleve  transcendent  (see,  for  example. 
Miles  1978) ,  and  would  be  qualitatively  similar  to  tne 
linear  solution  in  terms  of  Airy  functions.  With  this  in 
mind,  we  study  the  development  of  tne  wavefield  from  tne 
leading  edge  of  a  pair  of  straight  caustics  caused  oy 
increasing  water  depth  in  section  6.7,  and  also  look  at  the 
situation  of  waves  reflected  by  a  distribution  of  following 
current  in  section  6.8. 

Finally,  in  section  6.9,  we  study  tne  combined 
refraction  and  diffraction  of  an  incident  wave  field  oy  tne 
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rip  current  velocity  distribution  given  by  Artnur  (1950) . 
Here,  the  nonlinearity  occurs  at  too  large  an  Ursell  number 
for  Stokes  theory  to  be  valid,  and  we  employ  instead  a 
dispersion  relation  based  on  the  empirical  findings  of 
Walker (1976) ,  as  mentioned  in  section  3.5. 


u. 


.  w- , 
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6.2.  Numerical  Approximation  of  the  Parabolic  Equation 

The  parabolic  equations  of  section  4.4  are  nonlinear 
and  in  general  have  variable  coefficients  due  to 
inhomogeneity  of  the  physical  domain  (in  terras  of  depth  and 
current  variations) .  Consequently,  a  numerical  treatment  is 
required  in  order  to  obtain  solutions  to  any  but  the 
simplest  problems.  The  parabolic  nature  of  the  equations 
leads  to  the  choice  of  a  solution  technique  based  on  the 
implicit,  forward-stepping  scheme  of  Crank  and  Nicolson, 
detailed  below.  In  analogy  to  heat  transfer  and  diffusion 
problems  which  are  also  parabolic  in  nature,  we  choose  the  x 
coordinate  as  the  forward-marching,  time-like  dimension, 
while  the  y-coordinate  then  serves  as  the  physical  spatial 
coordinate.  For  all  examples  studied  in  this  chapter,  we 
will  restrict  our  attention  to  a  rectilinear  coordinate 
system  as  shown  in  Figure  6.1.  The  solution  scheme  also 
requires  the  specification  of  lateral  boundary  conditions  on 
the  y-coordinate,  which  will  be  discussed  in  section  6.3,. 

For  waves  in  the  absence  of  currents,  we  will 
consider  both  the  lowest  order  approximation  (4.4.4), 
referred  to  hereafter  as  the  "Radder"  approximation,  and  the 
higher  order  approximation  (4.4.38),  referred  to  as  the 
"Booij"  approximation.  For  waves  in  the  presence  of 
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currents,  the  higher  order  approximation  is  quite  complex 
and  its  use  does  not  produce  significantly  different  results 
than  the  model  equation  (4.4.35)  in  the  examples  studied  in 
sections  6.8  and  6.9.  Consequently,  we  will  restrict  our 
attention  to  the  simpler  model  (4.4.35).  The  current- free 
equations  (4.4.4)  and  (4.4.38)  will  be  treated  separately 
from  the  wave-current  model  (4.4.35),  due  to  slight 
differences  in  finite-difference  approximations  of  the 
coefficients . 


Both  approximations  (4.4.4)  and  (4.4.38)  may  be 
represented  by  the  general  equation 


Ax  (*.-*)  A  +_L.fap)xA  +A(p\\> 

2  bp  fcp 


P,  A  (pA*  \  +  A  lA'zA  +  cu  w  A 

79 1 


(6.2.1) 


where 
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and  where  primes  (')  have  been  dropped  from  the  amplitude  A 
for  convenience.  Note  that  the  reference  phase  function  *•? 
is  still  given  by  ik6x.  The  Radder  approximation  is 
recovered  by  the  choice  P(  *0,  Pr*l/2,  while  the  choice 

Pj*l/4,  PJ^-3/4  gives  the  Booij  approximation. 

A  rectangular  grid  with  uniform  spacing  in  the  x  and 

y  directions  is  established  as  shown  in  Figure  6.1.  We 

define 

()aX  +  j  /  £  x  £  ft 

For  all  examples  studied,  the  depth  will  be  taken  as  being 

uniform  on  i*l,  and  the  wavenumber  on  the  first  grid  row 

will  be  chosen  as  k0.  Values  of  yj  are  defined  according  to 

■  Ij-')*1}  +13i-  j  1  4  w 

In  (6. 2. 5-6),  xr  and  yr  are  arbitrarily  chosen  to  allow  a 
shift  in  origin. 

The  parabolic  model  (6.2.1)  is  written  in  finite 
difference  form  using  the  implicit  Crank-Nicolson  scheme. 


The  coefficient  ^3  is  written  as 


>*X  C  0?p)J*  -t  0*p)j  J 


where  the  i  row  contains  known  values  of  A  and  the  i+1  row 
is  the  row  to  be  determined  by  the  implicit  step.  The 
second  derivative  term  (pAu)w.  is  approximated  according  to 


ip^\t '  (pu  i£>0&!  -Al  1 1 


Using  (6. 2. 7-8)/  the  parabolic  model  is  approximated  by  the 
finite-difference  equation: 


*  (5?$f  J  ^  J  •> 


-J^U 

Sp>? 
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Equation  (6.2.9)  may  be  arranged  so  as  to  solve  for  the  A  A<': 

based  on  the  known  values  of  A4.  Using  the  two-point 

lateral  boundary  conditions  discussed  in  section  6.3  leads 

to  a  tridiagonal  system  of  N  equations  in  the  N  unknown 

AT  ,  and  solutions  are  obtained  rapidly  using  the 

4 

double-elimination  scheme  described  by  Carnahan,  Luther  and 
Wilkes (1969) ,  with  their  numerical  algorithm  modified  to 
handle  complex  arithmetic.  Previous  to  this  step,  however, 
an  estimate  must  be  obtained  for  the  unknown  value  *A  j*'  in 
the  nonlinear  term  at  row  (i+1).  This  value  is  provided  by 
an  intermediate  explicit  step  given  by 


XT -a!  At  * 
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Note  that  in  the  intermediate  step  the  derivative  (kp)^  and 

-*s  A 

the  term  f j  are  off-center  with  respect  to  the  value  of  A, 
while  the  term  (pA^)^x  is  neglected. 

The  implicit  scheme(6 . 2. 9)  has  truncation  errors  of 
0  (  Ux)  ■  ,  (ay)*” )  .  While  it  is  generally  difficult  to  assess 
the  stability  of  nonlinear  difference  equations  with 
variable  coefficients,  it  is  known  that  the  linear,  constant 
coefficient  counterpart  of  (6.2.9)  is  unconditionally  stable 
with  respect  to  choices  of  Z  x  and  XL.  y.  Further,  the  scheme 
in  the  linear,  constant  coefficient  case  preserves  mean 
square  quantities  (here  related  to  the  wave  energy  /Aj  ~) 
identically,  subject  to  the  accuracy  of  the  lateral  boundary 
conditions. 

The  explicit- implicit  scheme  described  by  (6.2.9-10) 
was  found  to  be  stable  in  all  experiments  for  all  values  of 

l  .y,  i-H  I 

kh,  with  the  explicit  value  I  A j  ,  generally  falling  within 

I—  • 

10%  of  the  final  value  |Aj‘j  .  A  somewhat  different  scheme 
involving  a  two-step  implicit  iteration,  as  used  by  Kirby 
and  Dalrymple(1983b) ,  was  found  to  develop  an  unstable 


divergence  in  alternate  iterations  for  the  value  of  Aj  for 
values  of  kh  less  than  0.15  in  general.  Both  schemes 
exhibit  stable  behavior  for  ranges  of  parameters  where 
Stokes'  theory  is  valid. 


The  finite-difference  scheme  for  the  wave-current 
interaction  model  (4.4.35)  is  given  by 


+  JL  \tJ*> 
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n 

(f.2 .//) 


The  intermediate  value  A  j  is  provided  by  the  explicit  step 


3 


M 


This  scheme  is  used  in  the  examples  of  sections  6.8  and  6.9 


For  each  scheme,  linear  results  may  be  obtained  by 
neglecting  the  cubic  nonlinear  terms  and  eliminating  the 
explicit  steps  except  in  the  case  where  waves  are  breaking, 
where  the  explicit  step  is  used  to  obtain  the  forward  value 
“X  ^  used  in  estimating  the  dissipation  coefficient'  w.  This 
procedure  is  discussed  in  section  6.4. 

The  numerical  procedure  requires  the  specification 
of  a| ,  1< j<N  on  the  first  grid  row.  In  the  cases  studied 
below,  we  will  restrict  our  attention  to  monochromatic 
waves,  in  which  case  A*  may  be  specified  simply  by 


=  A0C 


I*  j  ‘N 


where  A0  is  the  initial  amplitude  and  J  is  the  angle  of 
incidence  of  the  wave  with  respect  to  the  x-direction.  In 
principle,  the  method  may  be  used  for  a  directional  spectra 
of  incident  waves  which  may  be  specified  by 


where  n  is  the  number  of  components  each  having  amplitude 
A^/  direction  ^  and  phase  shift  ,  with  the  restriction 
that  each  component  be  a  plane  wave  with  absolute  frequency 
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6.3.  Treatment  of  the  Lateral  Boundary  Conditions 

The  solution  of  the  numerical  problem  presented  in 
section  6.2  requires  the  specification  of  lateral  boundary 
conditions  on  the  y-boundaries  of  the  computational  grid. 
In  the  examples  studied  in  this  chapter,  two  types  of  simple 
boundaries  are  included;  completely  reflective,  and 
completely  transmitting .  The  completely  reflective  barrier, 
representative  of  a  vertical  wall,  is  considered  first.  For 
a  plane  wall  oriented  in  the  x-direction,  the  condition  of 
no  velocity  component  normal  to  the  wall  leads  to  the 
condition 

Aj  =o 

*»  *  'Jur 

where  y  is  the  position  of  the  wall.  For  all  cases 
considered  below,  the  wall  will  be  centered  midway  between 
two  adjacent  grid  rows,  as  shown  in  Figure  6.2a.  The  finite 
difference  form  of  (6.3.1)  centered  on  the  wall  is  given  by 

A^+,  »  Aj  (c.iz) 

where  the  wall  is  located  at  yw  *(j+l/2)  —  y.  Taylor  series 
expansion  of  (6.3.2)  about  y^  leads  to  the  result 
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indicating  second  order  (  in  Ay  )  accuracy  for  the  boundary 
condition.  The  centering  of  the  wall  between  two  adjacent 
grid  rows  is  accomodated  in  the  schemes  of  section  6.2  by 
generating  auxiliary  grid  rows  to  handle  the  required 
overlapping  of  wave  data  on  the  two  rows  adjacent  to  the 
wall . 

Open,  transmitting  boundaries  may  occur  on  either 
the  up-  or  downwave  side  of  the  computational  grid.  For  the 
downwave  boundary,  the  boundary  condition  must  allow  waves 
to  radiate  freely  from  the  numerical  grid  with  as  little 
reflection  as  possible,  while,  for  the  upwave  boundary,  the 
condition  must  allow  the  entry  of  a  wave  to  replace  the 
waves  propagating  away  from  the  boundary  into  the  grid 
interior,  thus  avoiding  the  creation  of  a  shadow  adjacent  to 
the  grid  boundary. 

The  application  of  an  upwave  boundary  is  in  general 
fraught  with  difficulty,  since  no  definite  information 
beyond  pure  assumption  exists  to  indicate  the  exact  nature 
of  waves  approaching  the  boundary  from  the  exterior  of  the 
grid.  For  the  application  to  the  study  of  a  shore  attached 
breakwater  on  a  plane  beach  in  section  6.6,  the  lateral 
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boundaries  are  placed  far  from  the  obstacle,  and  we  may 
assume  as  a  first  step  that  the  wave  field  consists  entirely 
of  the  shoaling  and  refracting  incident  wave  and  the  wave 
scattered  from  the  breakwater.  The  incident  wave  has  the 
characteristics  of  uniform  amplitude  in  the  y-direction  and 
a  longshore  wavenumber  component  m  which  is  constant  as  a 
result  of  Snell's  law.  If  we  presume  that  the  boundaries 
are  far  enough  from  the  breakwater  such  that  no  scattered 
wave  reaches  the  boundary,  then  both  the  up-  and  downwave 
boundary  conditions  can  be  specified  exactly  by 

=  X* t  A  j  Cj  *  }  0.3.y) 

Once  again,  assuming  that  the  physical  grid  boundary  is 
centered  between  adjacent  grid  rows,  as  in  Figure  6.2b, 
(6.3.4)  may  be  rewritten  in  finite  difference  form  as 


where  m  is  presumed  to  be  known  precisely.  (6.3.5)  is 
rearranged  to  read 

Aj^  ( I  -  )  *  Aj  +  ft. 3.^) 


Introducing  Taylor  expansions  of  A  centered  at  y=(j+l/2)Ay 


leads  to  the  result 


4^  ~  XokA  +  /i  =  si*.  A  +  0(*j) 

lH  (f.3.7) 

The  errors  involved  in  (6.3.5)  are  thus  seen  to  be 
comparable  in  order  to  the  error  encountered  in  (6.3.2); 
however,  the  larger  coefficient  and  lower  order  derivative 
in  (6.3.7)  as  compared  to  (6.3.2)  indicated  that  larger 
numerical  errors  are  likely  to  be  associated  with  the 
approximate  radiation  condition  (6.3.5). 

For  cases  where  the  value  of  m  cannot  be  determined 
based  on  Snell's  law,  or  where  the  scattered  wave  represents 
a  significant  contribution  to  the  total  wave  field,  the 
value  of  m  must  be  estimated  on  a  numerical  basis.  For  the 
case  of  a  wavefield  consisting  of  more  than  one  component, 
we  assume  that  the  wave  may  be  estimated  locally  at  the 
boundary  as  a  plane  wave  of  unknown  direction  and  slowly 
varying  amplitude  A  in  the  y-direction.  Then  (6.3.4)  may 
be  used  again.  Assuming  that  the  wave  field  changes  slowly 
in  the  x  and  y  directions,  m  at  grid  row  x*i^x  may  be 
estimated  by  the  results  at  grid  row  x*(i-l)Ax  by 

*1*  =  -si  A^~' / A'" 
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Then,  applying 

the  centered 

finite  difference 

scheme  at 

x* ( i-1 ) A  x ,  y= 

( j+1/2 )  ^y 

yields  the  estimate 

for  the 

numerical  value 

m 

'""JJ 

•  -u  ( 

A.\+< 

-  *r ) 

(u.i) 

^  \ 

A+- 
n  j+( 

This  estimate  introduces  a  further  numerical  error  which  may 
J  be  estimated  by  inspecting  the  case  where  m  is  known  and  Ja| 

is  uniform  in  the  x  and  y  directions.  Employing  Taylor 

i 

i  series  expansions  about  the  center  point  yMj+l/ZJ^y  leads 

to  the  result 

I 

[  'hL*  *  ->  (  Au,/A  +  A?) 

I  (l  +  Av> ) 

I 

I 

1  Then,  assuming  that  the  relation  (6.3.4)  is  exact,  (6.3.10) 

! 

|  reduces  to 

!  ni  -  •*.  ( i  +  )  fa") 

j  w  it 

I 

I 

which  again  has  an  0(  ^y")  error  due  to  the  differencing 
scheme.  Substituting  (6.3.10)  into  the  error  estimate  for 
the  boundary  condition  (6.3.7)  leads  to  the  revised  formula 
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indicating  that  the  approximation  is  actually 
self-correcting  at  Retention  of  higher  powers  in 
Ay  in  each  expansion  shows  that  the  discretization  errors 
actually  cancel  to  any  order  in  Ay,  rendering  the  boundary 
condition  exact  in  regions  where  Snell's  law  holds  and 
amplitude  is  uniform  in  the  y  direction.  However  errors  may 
be  expected  to  occur  in  short  crested  wave  fields ,  where  A 
may  have  large  variations  over  several  grid  spacings. 

Results  of  several  numerical  experiments  are 
presented  here  in  order  to  indicate  the  errors  associated 
with  the  application  of  the  radiating  boundary  condition  to 
both  up-  and  downwave  open  boundaries.  Tests  were  conducted 
using  a  single  wavelength  and  period  for  waves  propagating 
over  a  flat  bottom.  Water  depth,  wavelength  and  period  were 
given  by 


h  *  10m. 

L  *  100m.  ;  k  ■  0.0628m 
T  »  10.726s. 


The  formulation  (6.3.6)  was  tested  using  both  the 
exact  value  m«ksin  s  and  the  numerical  value  m^  ,  given  by 


(6.3.9),  for  angles  c  *  10,  20,  30,  and  40  ,  where  o  is 
the  angle  between  the  direction  of  wave  propagation  and  the 
x>axis,  and  for  grid  spacings  (  ix,  4y)/L  *  0.1,  0.2,  0.3, 
0.4,  and  0.5.  Typical  errors  associated  with  the 
discretization  of  the  spatial  domain  occur  in  the  form  of 
partial  reflection  of  the  incident  wave  at  the  downwave 
boundary,  and  decay  of  wave  energy  at  the  upwave  boundary, 
again  due  to  reflection  of  waves  trying  to  enter  the  domain. 
For  tests  using  the  exact  value  of  m  and  the  monochromatic 
wavefield,  the  reflection  coefficient  at  the  downwave 
boundary  and  the  loss  of  energy  at  the  upwave  boundary  are 
plotted  as  functions  of.  wave  angle  and  grid  spacing  in 
figures  6.3  and  6.4,  respectively.  The  loss  of  energy  is 
represented  as  a  percent  error  between  the  minimum  energy 
occuring  within  ten  wavelengths  of  the  starting  grid  row  and 
the  starting  value.  Calculations  were  performed  using  the 
difference  scheme  (6.2.9)  in  linearised  form  with  Pj  *0. 
Results  indicate  that  the  numerical  errors  are  relatively 
small  for  all  angles  of  incidence  using  small  grid  spacing, 
and  increase  rapidly  with  the  loss  of  detail  in 
representation  of  the  wave  at  higher  grid  spacings. 

Results  using  the  numerically  exact  condition 
employing  the  calculated  value  of  m  verified  the  condition 
as  being  an  exact  representation  of  the  boundary  conditions 


for  plane  waves.  No  errors  on  either  the  upwave  or  downwave 
boundary  were  detected  within  the  limits  of  the  numerical 
accuracy  of  the  calculations  for  any  grid  spacing  or  angle 
of  incidence.  The  exactness  of  the  numerical  boundary 
condition  coupled  with  energy  conserving  properties  of  the 
Crank-Nicolson  scheme  allowed  the  use  of  grid  spacings  in 
excess  of  a  wave  length  with  no  loss  in  energy  in  the 
computational  domain,  even  though  the  individual  waves 
cannot  be  resolved  at  this  numerical  scale.  A 
representative  wavefield  with  4x/L=0.2,  6 *45*  is  shown  in 
Figure  6.5.  The  exactness  of  the  boundary  condition 
employing  a  numerical  value  m£  was  further  verified  by 
testing  the  case  of  waves  shoaling  on  a  plane  beach,  where 
Snell's  law  is  applicable,  and  again  no  errors  associated 
with  the  open  boundary  conditions  were  detected. 


As  a  final  test  of  the  applicability  of  the  open 
boundary  scheme  to  a  general  wavefield,  several  runs  were 
performed  for  the  case  of  two  intersecting  wave  trains,  both 
travelling  at  an  angle  to  the  grid  direction  x.  Both  waves 
were  taken  to  have  the  same  wavelength  as  the  wave  in  the 
single  wave  test,  and  one  wave  was  oriented  at  an  angle  of 
10 J  for  each  test  run,  with  the  angle  of  the  second  wave 
being  varied  from  15'’  to  40c  .  Grid  spacings  were  varied  as 
in  the  single  wave  tests.  A  representative  computed  wave 


field  is  shown  in  Figure  6.6  for  the  case  of  Ay/L  *  0.1  and 
9  *  20°  for  wave  2. 

The  distortion  of  the  wavefield  adjacent  to  the  up- 
and  downwave  boundaries  (left  and  right  boundaries, 
respectively,  in  the  figure)  is  a  result  of  the  assumption 
that  the  wavefield  could  be  characterised  by  a  plane  wave 
near  the  boundaries.  On  the  upwave  boundary,  the  boundary 
condition  assesses  the  characteristics  (amplitude  and 
composite  direction  of  propagation)  of  the  wave  at  the 

e 

incident  (x*0)  boundary  and  tnen  lets  a  plane  wave  with 
these  characteristics  propagate  into  the  domain,  creating  a 
plane  wave  in  the  geometric  region  shadowed  by  the  starting 
point  (x«0)  of  the  side  boundary.  On  the  downwave  boundary, 
the  short-  crested  pattern  propagates  out  of  the  domain  with 
little  difficulty  until  a  nodal  line  approaches  the  boundary 
(about  1/3  of  the  total  boundary  length  from  the  top) .  At 
this  point,  the  gradient  of  surface  elevation  normal  to  the 
boundary  is  dominated  by  the  short  crested  pattern  rather 
than  by  the  wavenumber  component  normal  to  the  boundary,  and 
the  assumptions  used  in  deriving  the  boundary  condition 
again  break  down.  The  model  responds  by  reflecting  a 
spurious  wave  back  into  the  grid  which  can  be  seen  through 
the  disruption  of  the  short-crested  wave  pattern  at  a 
considerable  distance  into  the  wavefield  at  the  terminating 
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boundary.  This  effect  can  become  severe  for  large  angles  of 
incidence  with  respect  to  the  x-axis. 

Although  the  basic  computational  scheme  developed  in 
section  6.2  is  clearly  capable  of  modelling  the  propagation 
of  a  directional  spectrum  of  waves,  it  is  clear  that  more 
development  of  open  boundary  conditions  must  be  pursued  in 
order  to  develop  a  model  of  a  fully  general  nature. 


6.4.  A  Model  for  Breaking  Waves  m  the  Surfzone 


An  advantage  of  the  parabolic  method  outlined  above 
over  solution  techniques  for  elliptic  and  hyperbolic 
equations  is  that  no  downwave  boundary  condition  is  needed 
for  the  solution  of  the  initial  boundary  value  problem. 
However,  in  applications  of  wave  models  to  coastal  areas, 
the  behavior  of  waves  in  the  vicinity  of  a  physical  downwave 
boundary  consisting  of  the  actual  coastline  is  of  primary 
importance  to  the  prediction  of  known  physical  effects  such 
as  the  wave  induced  setup  and  longshore  currents. 


Wave  breaking  in  the  surfzone  is  a  complex,  highly 
nonlinear  phenomenon.  It  is  obvious  that  the  model 
described  in  this  study,  which  is  limited  to  the 
representation  of  weak  nonlinearities,  is  basically 
incapable  of  representing  the  underlying  physics  of  the 
breaking  process.  However,  some  progress  can  be  made  by 
shifting  our  view  of  the  model  from  its  physical  basis  to 
its  use  as  a  predictive  tool. 


The  forces  leading  to  the  generation  and 
maintainance  of  setup  and  wave-^iftduc^d  currents  depend  on  a 
physical  balance  between  gradients  of  excess  momentum 
fluxes,  pressure  forces  due  to  changes  in  mean  surface 


the  wave  field.  Thus,  as  a  lowest  approximation  of  the 
overall  physics,  it  suffices  that  the  wave  model  be  able  to 
predict  ti:e  local  wave  amplitude  in  the  breaking  zone  with 

e 

some  degree  of  reliability.  The  simplest  model  of  wave 
decay  in  the  surfzone  is  based  on  the  assumption  that  the 
ratio  of  wave  height  to  local  water  depth  has  the  same  value 
everywhere  in  the  surfzone  as  at  the  breaker  line.  This 
assumption  has  been  used  extensively  in  the  literature,  from 
the  earliest  predictions  of  setup  (  Longuet-  Higgins  and 
Stewart  (1963),  Bowen,  Inman  and  Simmons ( 1968 ) )  and 
longshore  currents  (  Bowen(1968),  Longuet-Higg ins ( 1970) )  up 
to  the  latest  applications  of  numerical  refraction  schemes 
to  the  study  of  wave-induced  circulation  over  arbitrary 
bottoms  (  Ebersole  and  Dalrymple (1980) ,  Wu(1983)  ). 
However,  it  has  long  been  known  that  breaking  waves, 
especially  of  the  plunging  type,  do  not  follow  so  simple  a 
rule.  Extensive  model  tests  of  normally  incident  wave 
trains  breaking  on  laboratory  beaches  have  shown  that  the 
pattern  of  wave  height  decay  across  the  surfzone  is  strongly 
a  function  of  the  beach  slope.  Representative  measurements 
of  Horikawa  and  Kuo  (1966)  are  shown  for  example  in  Figure 
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The  purpose  of  the  present  section  is  to  relate  an 
empirical  model  of  surfzone  wave  energy  decay  to  the 
coefficient  w  of  the  wave  model,  and  to  detail  the 
application  of  the  model  to  the  prediction  of  wave  height  in 
the  surfzone.  The  empirical  model  is  taken  from  the  work  of 
Dally(1980) . 


Dally  proposed  that  the  decay  of  energy  flux  with 
distance  in  the  surfzone  ba  given  by  the  relation 

(Ec>)x  -  -&  (ec^  -  )  (i.1.1) 


where  h  is  the  local  water  depth  and  K  is  a  constant  to  be 
determined,  and  is  related  to  the  rate  of  energy  decay.  The 


quantity  (ECg)s  is  a  "stable"  energy  flux  for  a  broken  wave 
in  water  of  depth  h.  The  stable  energy  flux  may  be  related 


to  the  height  obtained  by  a  wave  propagating  over  a  flat 


bottom  after  the  cessation  of  breaking.  Data  obtained  by 


Horikawa  and  Kuo  (1966)  suggest  that  waves  typically  reach  a 
height  of  0.4h  at  the  cessation  of  breaking  on  a  flat 
bottom,  while  an  asymptotic  value  of  0.5h  is  approached  for 


waves  breaking  on  a  plane  slope.  In  the  following 
derivation,  we  denote  (H/h)s  *  Z  ,  where  c  is  a  free 
parameter . 
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Equation  (6.4.1)  may  be  related  to  the  wave  energy 
equation  (E.ll)  after  neglecting  currents  and  assuming  a 
time-steady  wave  field.  For  one-dimension,  (E.ll)  becomes 

(EC%)X  -  -wC 

Noting  that  C^g*  C<j  ,  w  may  be  written  as 

(<•«> 

where  H  is  the  local  wave  height,  given  by  H=2,Ai.  Since  H, 
and  hence  w,  is  not  known  at  the  location  of  the  forward 
step  in  the  implicit  schemes  of  section  6.2,  an  iterative 
scheme  must  be  employed  to  obtain  an  estimate  of  w  at  the 
forward  position  before  the  final  application  of  the 
Crank-Nicolson  scheme. 

For  a  plane  slope  and  neglecting  the  effect  of  setup 
(which  is  not  calculated  in  the  wave  model) ,  a  simple 
analytic  solution  of  (6.4.1)  may  be  obtained  and  used  for 
comparison  with  the  numerical  model.  Letting  "J  3  EC^ 
denote  the  energy  flux  and  writing  h(x)  as  h»-s(x-x0) ,  where 
x0  is  the  position  of  the  shoreline,  (6.4.1)  may  be  written 


where  -i  »  K/s.  Assuming  shallow  water  conditions,  (EC^  )$ 
may  be  written  as 


(ecj )s  * 


ft.tr) 


Equation  (6.4.4)  then  has  the  general  solution 

3-*c^  * 


where  g^/8.  A  special  solution  is  needed  for  ci  •  5/2. 
The  value  of  c  may  then  be  determined  by  the  value  of  at 
the  breaker  line,  where 

k  *  hy  ;  H  *  »  "ft  )  X  * 


and  ~K  is  the  ratio  of  breaking  wave  height  to  still  water 
depth  at  the  breaker  line.  is  given  by 


3b  =  ^V/l 


which  leads  to 


c  •/**'[! -(***.)&)  b.n) 

,  z 

After  rewriting  cT  as  j  *  |J  H  h-“,  (6.4.6)  may  be  rewritten 


in  dimensionless  form  as 


where 


o<-4 


«.*4) 


For  the  case  C%  *  5/2,  the  method  of  variation  of 
parameters  and  application  of  the  boundary  condition  at  the 
breaker  line  yields  the  result 


Based  on  a  comparison  of  the  laboratory  data  of 
Horikawa  and  Kuo (1966),  Dally  chose  the  value  K*0.17.  The 
special  case  *  5/2  then  corresponds  to  a  beach  slope 
s*0.068.  Results  for  a  range  of  values  of  1  <  !  <  10, 
corresponding  to  the  range  of  beach  slopes  0.17  >  s  >  0.017, 
are  given  in  Figure  6.8  for  u  »  0.4,  “^  *  0.78.  The  lines 


i5» 


labelled  1  and  2  correspond  to  the  constant  decay 
H*"Kh»0.78h  and  to  the  stable  wave  height  Hs=^h  *  0.4h, 
respectively. 


The  results  (6.4.8-10)  provide  a  check  for 
determining  the  accuracy  of  iterative  schemes  using  the  wave 
damping  term  (6.4.3).  Noting  that  Hs  »  h  and  H*2jAj, 

(6.4.3)  may  be  rewritten  as 

v  *  kc.  r  i  -  yv  i  (m-u) 

h  TiTp- i 

We  use  the  linearized  form  of  the  parabolic  equation 

(4.4.4) ,  which  can  be  written  as 

Ax  -  A.  A  ***  ,*1—  A  +  yUT  A  *  O  Cb-  */.  /2.) 

where 

/o-  *  X  ,  Xr  l  -  jCJl  1  Ci-v.13) 

7  2fi  L  V  Ml1  J 

(6.4.12)  is  written  in  finite  difference  form  as 
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where 


A-+! 


*X(h*”)X  1 

H  l A** 1 1  J 


ft.V./r) 


The  intermediate  value  '~+!  is  determined  by  a  similar  step 


A*”- A'  4(k*QX*"  *  Ck'-Wl 


AX 


*i  L  C, 


x"+'  J-'1 


%  ch 


+  X  [jur*1  A*  O 

If)  ah) 


where 


A# 

A 


v^ZT  i 

via*  J 


Several  cases  were  run  for  waves  starting  in  a  depth 
of  2m  and  propagating  towards  shore  over  a  plane  slope.  The 
program  checks  at  each  step  that  the  wave  height  has  not 
exceeded  the  breaking  criterion.  When  H  becomes  greater 
than  *>.  h,  the  program  begins  calculating  values  of  the 
damping  coefficient  (6.4.15).  Breaking  continues  until^r 
falls  to  a  value  of  zero,  which  does  not  occur  on  the  plane 
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beaches  studied  here  but  would  be  expected  to  occur  readily 
for  waves  propagating  over  uneven  topography.  For  each 
case,  the  wave  height  was  assumed  to  be  1.0m  at  a  depth  of 
2m,  and  wave  period  was  assumed  to  be  5  seconds.  Values  of 
1,  3,  and  10  were  tested  using  various  computational 

grid  spaces  «x.  Results  for  ■*.*!  and  u  x-0 . 2m  and  1.0m  are 
shown  in  Figure  6.9,  3=3  and  tlx*  1.0,  2.0,  and  5.0m  in 

Figure  6.10,  and  oi  =10  and  <ix*  2.0,  5.0,  and  10.0m  in 

Figure  6.11.  The  exact  solution  (6.4.8)  is  included  in  each 
figure  for  comparison,  with  h^  being  taken  as  being  the 
average  of  the  depth  at  the  last  grid  point  before  breaking 
and  at  the  first  grid  point  after  breaking.  In  each  case, 
the  numerical  results  provide  an  adequate  representation  of 
the  exact  solution.  In  Figure  6.12,  the  entire  process  of 
shoaling  from  h=2m  up  to  the  break  point  and  the  subsequent 
decay  in  the  surfzone  is  shown  for  the  case  *10  and 
^x*5m,  with  the  solution  (6.4.3)  included  for  comparison. 

The  breaking  wave  model  has  been  incorporated  into 
each  of  the  numerical  schemes  discussed  in  section  6.2  by 
adapting  it  to  the  explicit- implicit  iteration  scheme,  and 
is  used  in  the  remainder  of  the  examples  in  Chapter  6  where 
appl icable. 

The  model  described  in  this  section  is  empirical  in 


r/rjO 
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nature  and  is  based  on  a  two-parameter  best  fit  to  existing 
laboratory  data.  An  alternate  procedure  based  on  estimating 
the  rate  of  energy  decay  in  a  turbulent  bore  has  been 
suggested  by  Divoky,  LeMehaute  and  Lin (1970)  and  has  been 
revised  to  include  the  case  of  periodic  waves  by 
Batt j es ( 1978 ) .  While  these  methods  attempt  a  more  thorough 
approach  to  the  physics  underlying  the  breaking  process, 
they  do  not  differ  with  the  model  described  here  in  terms  of 
results.  Indeed,  Divoky,  LeMehaute  and  Lin  based  their 
conclusions  on  the  relative  validity  of  their  model  on 
comparison  with  the  same  set  of  laboratory  data  used  to 
calibrate  the  present  model. 
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6.5.  Combined  Refraction  and  Diffraction  of  Waves  by  a 
Submerged  Shoal 


The  numerical  approximations  developed  in  sections 
6. 2-6. 4  lead  to  a  parabolic  equation  method  which  may  be 
used  to  model  wave  propagation  towards  an  open  coastline  up 
to  and  including  the  surfzone.  Various  sets  of  experimental 
data  exist  which  may  be  used  to  test  the  predictions  of  the 
nonlinear  model.  Recently,  Berkhoff,  Booy  and  Radder  (1982) 
have  presented  data  obtained  in  the  vicinity  of  a  focus  and 
cusped  caustic  created  by  an  elliptic  shoal  resting  on  a 
plane  slope,  and  have  compared  the  data  to  the  predictions 
of  three  linear  wave  models:  a  refraction  scheme  involving 
averaging  over  bundles  of  adjacent  rays  (Bouws  and  Battjes, 
1982),  a  parabolic  equation  model  for  the  scattered  incident 
wave,  and  an  elliptic  model  for  the  entire  wave  field. 
While  the  results  of  each  computational  model  differ  in 
particulars,  an  important  conclusion  may  be  drawn  with 
respect  to  comparison  between  linear  wave  models  and  the 
data  set.  Linear  models  uniformly  tend  to  overpredict 
maximum  wave  amplitudes  in  the  vicinity  of  focussed  waves, 
where  wave  steepness  may  become  large  and  nonlinear  effects 
become  important.  This  tendency  is  clear  in  the  comparison 
between  Berkhoff,  Booy  and  Radder' s  parabolic  model  results 
and  data.  The  conclusion  to  be  drawn  from  a  comparison  of 
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data  and  their  elliptic  model  is  less  clear.  The  elliptic 
model  tends  to  predict  lower  amplitudes  in  the  region  of 
focussed  waves  than  the  parabolic  model;  however,  the  model 
results  are  significantly  contaminated  with  waves  reflected 
from  the  downwave  boundary,  as  evidenced  by  the  rapid 
oscillations  in  amplitude  along  x-direction  transects  6-8  in 
their  results.  The  effect  of  these  reflected  components  on 
the  structure  and  development  of  the  the  focus  is  unknown, 
and  it  is  unclear  that  the  predictions  of  the  elliptic  model 
are  any  better  than  those  of  the  parabolic  model  in 
comparison  to  the  experimental  data. 


The  intent  of  the  present  section  is  two-fold. 
First,  we  wish  to  verify  that  a  parabolic  equation  for 
weakly-nonlinear  wave  motion  is  capable  of  predicting 
accurate  results  for  waves  in  the  vicinity  of  a  cusped 
caustic,  by  comparing  the  model  predictions  to  data. 
Secondly,  by  contrasting  the  difference  in  predictions  of 
linear  and  nonlinear  models,  we  wish  to  show  that  the 
difference  between  previous  linear  model  results  and  data  is 
largely  due  to  the  neglect  of  nonlinearity  rather  than  to 
any  inherent  inaccuracy  in  the  modelling  techniques.  We 
will  concentrate  on  the  data  set  of  Berkhoff,  Booy  and 
Radder  in  performing  this  comparison.  This  data  set  is  well 
suited  to  the  present  purpose  due  to  the  great  detail  in 


which  data  on  wave  amplitude  was  obtained  over  the  entire 
vicinity  of  a  refractive  focus,  where  diffractive  and 
nonlinear  effects  become  significant.  We  will  use  the 
"Bladder"  approximation  (4.4.4)  in  its  numerical  form 
(6.2.9,10)  in  this  section.  Computation  is  halted  before 
the  breaker  line  is  reached;  consequently  we  neglect  the 
effect  of  dissipation  due  to  breaking. 


Experimental  topography  and  computational  domain 


Details  of  the  experimental  procedure  and  setup  may 
be  obtained  from  Berkhoff (1982)  or  Berkhoff,  Booy  and 
Radder (1982) ,  which  include  a  photograph  of  the  experimental 
wave  field  and  a  plot  of  wave  rays  in  the  refraction 
approximation;  we  summarize  the  important  points  here.  The 
experimental  topography  consists  of  an  elliptic  shoal 
resting  on  a  plane  sloping  bottom  with  a  slope  of  1:50.  The 
plane  slope  rises  from  a  region  of  constant  depth  h=0.45m, 
and  the  entire  slope  is  turned  at  an  angle  of  20~  to  a 
straight  wave  paddle. 


The  bottom  contours  and  chosen  computational  domain 
are  shown  in  Figure  6.13,  where  the  computational  domain  is 
represented  by  the  dashed  line  bounding  the  contours. 
Experimental  data  is  given  for  the  sections  labelled  1-8. 


Computational  coordinates  are  established  with  the  origin  at 
the  upper  left  corner  of  the  domain,  which  is  then  given  by 
X^  21, vT 'm.  .  0<l  tj.  ^  ZO.O'm. 


O^x  (a,  £  ZO.O'm. 


The  offshore  boundary  of  the  domain  is  chosen  so  that  water 
depth  is  constant  along  x=0.  The  initial  condition  for  the 
wave  then  corresponds  to  the  uniform  wave  train  generated  at 


the  wave  paddle;  we  set 


A(x~o j  *  Ae 


where  Aft*0. 0232m  is  the  amplitude  of  the  incident  wave.  The 
wave  period  T=1  sec.  Slope-oriented  coordinates  {x',y‘}  are 
established  which  are  related  to  the  computational 


coordinates  {x,y}  according  to 


X*  s  (  X  -  /0.x')  op.  2J>*  -  to)  s i*>  2o* 
^  a  [v.  -  iox)  S:i>  Zo*  +  (oy-  io>:oj2o* 


(CX  3a 


((>.5, lie 


The  origin  {x’,y'}  =  {0,0}  corresponds  to  the  center  of  the 


shoal.  The  slope  is  described  by 


0.  HS  *«. 


0.hc*l  -  0.6-l  (x  12  +  x' )  *«.  .  x'h  fax  ft 


The  boundary  of  the  elliptic  shoal  is  given  by 


&.«■) 


and  the  depth  in  the  shoal  region  is  modified  according  to 

h  =  Ku?t  -  [°s  (I  -faf-terff  -  *3J  &«) 

resulting  in  a  depth  h(x’,y‘*0)  =  0.1332m. 

The  lateral  boundaries  at  y=0,20m  are  open,  and 
transmitting  conditions  are  specified  according  to  (6.3.5), 
where  we  have  assumed  that  |A,  varies  slowly  in  the  y 
direction  at  the  boundaries.  The  wavenumber  component  m  at 
each  grid  row  is  evaluated  numerically  according  to  (6.3.9). 

Comparison  of  Computational  Results  and  Experimental 
Measurements 

The  nonlinear  equation  is  valid  under  the  same 
conditions  as  the  Stokes'  theory,  with  the  principle 


condition  being  given  by 


Uh  =  pAj  ]/(kh)Z  <  O0) 


&.T.7) 


where  Or  is  the  Ursell  number.  For  conditions  where  0^. 
approaches  0(1),  wave  models  based  on  the  Boussinesq 
equations,  such  as  that  of  Abbott,  Petersen  and 
Skovgaard (1978)  become  more  appropriate.  Also,  the 
derivation  of  (4.4.4)  technically  requires  that  r7j,h/(k  [a,:  ) 
be  0(a),  although  this  has  not  been  found  to  be  a  necessary 

a 

restriction  in  applications  of  the  linear  model. 


Laboratory  experiments  designed  to  test  the 
predictions  of  linear  wave  models  typically  satisfy  the 
condition  of  small  0^.  in  the  vicinity  of  the  wavemaker . 
However,  as  waves  propagate  into  shallower  water  and  are 
focused  by  bottom  irregularities,  0^.  may  increase  rapidly, 
and  care  must  be  taken  in  verifying  that  Stokes'  wave  theory 
is  valid  for  the  problem  being  considered. 


Preliminary  studies  of  the  relevant  parameters  in 
the  wave  field  were  performed  by  determining  Ur  at  several 
locations  in  the  experimental  wave  field.  Values  of  Ur  at 
the  wavemaker,  shoal  crest,  and  point  of  maximum  amplitude 
in  the  focused  region  are  given  by  0 ^  *  .014,  .290,  and 


.213,  respectively.  Stokes'  theory  should  thus  be  valid 
throughout  the  region  of  principle  interest. 

Data  from  the  laboratory  experiment  of  Berkhoff, 
Booy  and  Radder  (1982)  is  available  for  the  labelled 
sections  1-8  indicated  in  Figure  6.13.  The  computational 
domain  was  discretized  into  square  grids  (  *1  x  =  *ly  *  grid 
spacing) ,  and  the  grid  scheme  was  established  so  that  grid 
rows  coincided  with  the  measurement  transects.  Grid  size 
was  decreased  until  the  point  was  reached  where  further 
reduction  did  not  affect  model  predictions  significantly. 
The  final  numerical  calculations  were  performed  using  a 
spacing  of  A  x  *  0.25m,  corresponding  to  a  grid  of  87X81 
rows . 

Contours  of  normalized  wave  amplitude  are  shown  in 
relation  to  the  bottom  contours  for  the  linear  wave  field  in 
Figure  6.14  and  for  the  nonlinear  wave  field  in  Figure  6.15. 
The  nonlinear  results  indicate  a  broadening  of  the  central 
focused  region  and  a  decrease  in  the  maximum  amplitude  in 
the  focus.  The  nonlinear  contours  are  in  better  agreement 
with  the  experimental  contours  shown  in  Figure  2  of 
Berkhoff,  Booy  and  Radder ( 1982 ) . 

Results  of  the  nonlinear  and  linear  models  are  shown 


$ 


On  section  1,  focussing  effects  have  only  begun  to 
be  apparent,  and  only  a  slight  difference  exists  between  the 
linear  and  nonlinear  models,  which  both  agree  well  with  the 
data.  Sections  2  and  3  describe  the  region  of  the 
development  of  the  cusped  caustic.  For  both  sections,  data 
typically  falls  between  the  predictions  of  the  two  models  in 
the  region  of  maximum  amplitude,  with  the  nonlinear  model 
underpredicting  the  maximum  amplitude  by  about  8%  on  section 
3.  However,  on  sections  4  and  5,  where  the  wave  has  passed 
through  the  cusped  caustic,  the  nonlinear  model  predictions 
and  the  data  are  in  striking  agreement,  with  both  the  height 
and  width  of  the  central  focussed  region  and  the  size  and 
shape  of  side  lobes  in  the  diffraction  pattern  being  very 
well  predicted. 


Agreement  between  the  nonlinear  model  and  the  data 
is  also  verv  aood  alona  f.he  lonaitudinal  sections  6-8.  On 
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in  amplitude  in  the  region  of  the  cusp  (with  respect  to 
linear  model  predictions)  and  the  slower  decay  of  amplitude 
towards  the  shoreward  end  of  the  transect  (large  x)  .  On 
section  6,  the  nonlinear  model  quite  accurately  predicts  a 
sharp  drop  in  amplitude  which  is  apparent  in  the  data  but 
totally  absent  in  the  results  of  the  linear  model. 
Berkhoff,  Booy  and  Radder  interpreted  the  discrepancy 
between  the  linear  model  prediction  and  the  data  in  this 
region  as  being  the  result  of  the  proximity  of  an 
amphidromic  point  in  the  wave  phase  surface  and  the  linear 
model's  inability  to  resolve  this  point.  Here,  we  see  that 
the  discrepancy  is  entirely  explained  in  terms  of  nonlinear 
effects  and  the  resulting  change  in  the  geometry  of  the 
focused  region. 


Differences  between  the  linear  and  nonlinear  results 
are  less  striking  along  section  8;  The  nonlinear  model  does 
reproduce  the  relatively  high  amplitudes  in  the  region  6  to 
9  meters  downwave  of  the  shoal,  in  comparison  to  the  linear 
model . 


Taken  as  a  whole,  it  is  apparent  that  the  results  of 
the  nonlinear  model  exhibit  closer  agreement  with  the 
experimental  data  than  do  the  results  of  the  linear  model, 
with  the  only  obvious  discrepancy  occurring  in  the  region  of 


■  v »  v  v*.. 


initial  development  of  the  wave  focus.  It  is  possible  that 
the  discrepancy  in  this  initial  region  of  growing  amplitude, 
seen  in  sections  2  and  3,  is  due  in  part  to  a  failure  of  the 
approximate  equations  to  allow  waves  to  focus  rapidly  enough 
to  obtain  a  realistic  amplitude.  This  discrepancy  is  seen 
to  be  short-lived,  as  evidenced  by  the  remainder  of  the  data 
set. 


'‘•'V  V*,  7 A.  *>  V, 
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6.6.  The  Wavefield  around  a  Shore-attached  Breakwater 

As  a  second  test  of  the  parabolic  models,  we  next 
study  the  wavefield  in  the  vicinity  of  the  shore  attached 
breakwater  described  in  Figure  6.17.  An  extensive  set  of 
data  for  the  wavefield  in  the  shadow  zone  downwave  of  the 
breakwater  has  been  given  by  Hales (1980)  for  a  number  of 
wave  periods,  amplitudes  and  angles  of  incidence.  A  less 
detailed  set  of  data  is  available  for  the  reflection  zone  on 
the  upwave  side  of  the  breakwater  in  Pantazaras(1979) .  The 
given  experimental  arrangement  has  been  investigated  using 
several  computational  techniques  for  linear  waves,  including 
the  finite  element  method  for  the  elliptic  formulation 
(Houston',  1981)  and  a  parabolic  equation  method  based  on 
curvilinear,  ray-fitted  coordinates  (Tsay  and  Liu,  1982).  A 
closed  form  solution  in  the  linear,  mild-slope  approximation 
has  been  provided  by  Liu,  Lozano  and  Pantazaras (1979)  and 
has  been  compared  to  the  experimental  data  by  Liu(1982),  who 
found  qualitative  agreement  between  the  linear  theory  and 
experimental  results. 

As  a  first  test  of  the  parabolic  models,  we  restrict 
attention  to  linear  wave  theory  and  investigate  the 
difference  in  predictions  between  the  linearized  "Radder" 


and  "Booij"  approximations.  The  computational  grid  is  again 
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rectilinear,  with  the  x-direction  oriented  parallel  to  the 
breakwater.  Thus,  waves  incident  on  the  grid  may  be 
propagating  at  angles  of  up  to  30  to  the  x-direction.  Two 
cases  presented  by  Liu(1982)  are  studied,  and  results  are 
presented  for  a  y-direction  transect  5  ft.  shoreward  of  the 
breakwater  tip.  In  case  (a) ,  the  wave  period  T*1  sec.  and 

..  a 

tne  incident  wave  angle  oo*20  ,  while  for  case  (b) ,  T=1.5 
sec.  and  j^*30s.  Results  obtained  using  the  linearized  form 
of  (6.2.9)  with  Pj  *0  and  P^  *1/4  are  shown  in  Figures  6.13 
and  6.19,  respectively,  for  cases  (a)  and  (b) .  The  results 
of  Liu,  Lozano  and  Pantazaras (1979)  are  included  for 
comparison. 


Results  obtained  using  the  "Radder"  approximation 
(P, *0)  indicate  an  inability  to  model  the  diffraction  fringe 
on  either  the  up-  or  downwave  side  of  the  breakwater.  This 
result  is  particularly  apparent  in  the  region  adjacent  to 
the  shadow  boundary,  where  reflected  wave  activity  is 
minimal  or  absent  in  comparison  to  the  analytic  results. 
Similar  results  have  been  obtained  by  Hasimoto( 1982) ,  who 
used  the  original  model  of  Radder (1979) .  in  contrast,  the 
"Booij"  approximation,  with  P,  *1/4,  is  seen  to  give  a  better 
indication  of  the  wave  field  away  from  the  shadow  and 
reflection  zone  boundaries,  and  indeed  may  overestimate  the 
size  and  number  of  fringes  in  comparison  to  the  analytic 


Figure  6.19  Normalized  amplitude  perpendicular  to 
breakwater 

-  linear  waves/  *  1/4. 

-  Liu.  Lozano  and  Pantazaras  (1979) 


results.  The  overall  agreement  of  results  obtained  using 
*1/4  with  the  analytic  results  is  favorable  in  comparison 
to  the  results  obtained  by  Tsay  and  Liu(1982),  which  tend  to 
underestimate  the  diffraction  fringe  effects.  (This 

underestimation  may  be  the  result  of  interaction  with 
lateral  boundary  conditions  rather  than  limitations  of  their 
approximation.)  Both  models  are  seen  to  accurately  predict 
the  characteristics  of  the  reflection  and  shadow  zones  in 
the  near  field  of  the  breakwater  with  the  exception  of  the 
shift  of  the  first  amplitude  maximum  towards  the  breakwater 
seen  in  the  "Radder"  approximation. 

The  inability  of  the  linearized  form  of  the  "Radder" 
approximation  (Pf  *0)  to  correctly  model  amplitude 
modulations  at  a  distance  from  the  obstacle  disturbing  the 
wavefield  has  consequences  on  the  applicability  of  that 
model  to  the  study  of  nearshore  waves,  where  longshore 
amplitude  modulations  can  be  instrumental  in  establishing 
cellular  circulation  patterns  in  the  surfzone  and  nearshore 
region  (Dalrymple(1975) ,  Liu  and  Mei(1976)). 

The  contours  of  amplitude  and  the  breaker  line  are 
given  in  Figure  6.20  for  the  conditions  of  test  T8  of 
Pantazaras(1979) ,  assuming  a  breaker  index  of  0.78  and  using 
the  "Radder"  mode'  .  We  ***  :e  that  the  breaker  line  would 


Figure  6.20  Wavefield  in  vicinity  of  shore-attached 
breakwater.  Parameters  for  test  T8  of 
Pantazaras  (1979). 

T  »  1  sec.,  0Q  *  20°,  AQ  =  0.0320  m 

linear  model  results,  =  0. 

-  amplitude  contours  relative  to  AQ 

-  position  of  breaker  line  based  on 

H. /h,  *  0.78. 
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undergo  more  on-offshore  displacement  with  distance  upwave 
from  the  breakwater  in  the  "Booij”  model  due  to  the 
increased  amplitude  modulation  in  the  diffraction  fringes 
predicted  by  that  model.  However,  both  models  are  adequate 
for  predicting  wave  activity  in  the  immediate  vicinity  of 
the  breakwater. 

We  turn  now  to  a  consideration  of  nonlinear  effects 
in  the  wavefield.  In  particular,  we  wish  to  determine  if 
any  Mach-stem  effect  is  present  in  the  reflected  wavefield 
on  the  upwave  side  of  the  breakwater.  Nonlinear  effects  on 
the  downwave,  shadow  side  of  the  breakwater  would  be 
expected  to  be  confined  to  a  slight  increase  in  wave 
amplitude  in  the  shadow  region  due  to  increased  diffraction 
across  the  shadow  boundary;  see  Yue(1980),  section  4.3.  We 
choose  conditions  corresponding  to  tests  T5-T6  and  T8  of 
Pantazaras(1979)  in  order  to  perform  a  comparison  with  data. 
The  test  conditions  are  summarized  in  Table  6.1.  The 
comparison  is  performed  using  the  "Radder”  model  (P,  *0) ,  and 
we  compare  with  data  only  in  the  immediate  vicinity  of  the 
breakwater.  Tests  T5  and  T6  represent  identical  conditions 
of  wave  period  and  angle  of  incidence  and,  unfortunately, 
Pantazaras  presented  results  for  these  tests  in  averaged 
form,  since  he  was  only  interested  in  the  linear 
approximation.  Here,  we  will  compare  each  test  to  the 
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average  of  the  experimental  data.  Results  for  test  T8  ara 
shown  in  Figure  6.21.  Along  the  transect  at  x*2.5ft., 
nonlinear  effects  are  weak  and  there  is  no  apparent 
difference  between  the  two  solutions.  At  x»6ft.,  there  is 
some  evidence  of  Mach-stem  formation  in  the  nonlinear 
solution,  and  the  amplitude  adjacent  to  the  breakwater  is 
reduced  by  9%  in  comparison  to  the  linear  solution.  The 
experimental  data  is  somewhat  more  supportive  of  the  linear 
solution.  For  the  x*6ft.  transect,  the  (Jrsell  number  for 
the  wave  adjacent  to  the  breakwater  is  Up*0.648. 

Results  for  test  T5  are  presented  in  Figure  6.22. 
Again,  there  is  little  effect  due  to  nonlinearity  at  the 
x*2.5ft.  transect.  At  the  x**7.5ft.  transect,  the  tendency 
towards  the  formation  of  a  Mach  stem  is  apparent.  in  this 
case,  the  data  of  Pantazaras(1979)  is  more  supportive  of  the 
nonlinear  solution.  Here,  the  (Jrsell  number  for  the 
nonlinear  solution  at  the  breakwater  is  Up*1.99,  and  Stokes 
theory  is  of  questionable  validity  in  this  region.  However, 
the  large  value  of  Uj.  and  the  apparent  formation  of  a 
Mach>stem  indicate  the  importance  of  nonlinearity  in  this 
region.  Waves  in  this  region  would  be  expected  to  be  of 
cnoidal  form;  waves  of  this  class  also  are  known  to  produce 
Mach  stems  during  reflection  at  small  angles  of  incidence. 


Figure  6.21  Amplitude  contours  perpendicular  to 
shore-attached  breakwater.  Test  T8 
T  *  1  sec.,  8  *  20°,  AQ  *  0.0320  m. 

— —  linear;  -  nonlinear; 

•  Pantazaras  (1979)  data 


Figure  6.22  Amplitude  contours  perpendicular  to  shore- 
attached  breakwater.  Test  T5;  T  *  1.5  sec 
9  «  30°,  A  -  0.0229  m. 

— —  linear;  -  nonlinear; 

•  Pantazaras  (1979)  data. 


For  both  of  the  tests  T6  and  T7,  which  represent 
larger  amplitude  counterparts  of  tests  T5  and  T8 
respectively/  results  of  the  present  study  indicate  that 
wave  breaking  occurs  before  the  position  of  the  second 
transect  is  reached,  based  on  a  breaker  index  of  0.78.  In 
particular,  for  test  T6,  the  value  of  H/h  adjacent  to  the 
breakwater  from  the  linear  solution  without  breaking  is 
1.20.  Values  of  H/h  this  large  may  occur  for  waves  shoaling 
on  steep  beaches.  The  broken  wave  amplitude  resulting  from 
the  present  calculations  with  a  breaker  index  of  0.78  are 
presented  in  Figure  6.23.  We  note  that  the  Ur sell  number 
for  this  wave  condition  is  too  large  for  the  validity  of 
Stokes  theory. 

The  present  numerical  calculations  could  be  refined 
by  incorporating  a  variable  breaking  index  in  the  model. 
Various  formulae  exist  for  calculating  Hb/hb  based  on  the 
local  bottom  slope  and  parameters  of  the  wavefield.  This 
body  of  theory  and  empirical  results  has  been  reviewed 
recently  by  Mallard ( 1978) .  This  extension  to  the  model  is 
clearly  warranted  as  indicated  by  the  discrepancy  between 
position  of  the  breaker  line  calculated  here  and  estimated 


from  Pantazaras'  data. 
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Figure  6.23  Amplitude  contours  perpendicular  to 
shore-attached  breakwater.  Test  T6; 


T  *  1 . 5  sec. ,  0 


0.0344  m 


Values  based  on  breaking  index  <  *  0.78 

-  nonlinear;  -  linear 

•  Pantazaras  (1979)  data. 
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€.7.  Reflection  from  a  Caustic:  Topographic  Case. 

As  a  final  example  for  steady  waves  in  the  absence 
of  a  current,  we  study  the  incidence  of  a  plane  wave  in 
constant  water  depth  on  a  symmetric,  wedge-shaped 
depression,  with  sides  sloping  down  from  the  constant  depth 
region.  The  geometry  of  the  wedge  is  described  in  Figure 
6.24.  The  line  of  symmetry  is  taken  as  y»0,  and  the 
boundary  of  the  depression  is  given  by 

l ( i.7.1 ) 

where  £  is  the  wedge  half-angle,  and  where  the  tip  of  the 
wedge  is  located  at  x«0.  The  bottom  slopes  down  with  a 
slope  of  1:50  normal  to  the  wedge  boundary.  The  maximum 
depth  of  the  depression  is  given  by  2h# ,  where  h,  is  the 
depth  of  the  incident  wave  region.  The  incident  wave  is 
characterised  by  koho«1.0  in  the  constant  depth  region 
outside  the  wedge.  For  the  cases  considered  here,  a  caustic 
of  the  linear  wave  field  occurs  on  the  sloping  wedge 
boundary,  and  ,  in  the  far  field  (x  large) ,  the  wave  field 
in  the  vicinity  of  the  caustic  would  be  described  by  the 
Airy  function,  with  an  exponentially  decaying  amplitude  in 
the  geometric  shadow  over  the  trench  and  total  reflection  of 
the  incident  wave.  The  corresponding  result  for  weakly 
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nonlinear  waves  would  be  a  qualitatively  similar  description 
in  terms  of  the  second  Painleve  transcendent  (Miles  (1978)). 
However,  as  in  the  case  of  the  development  of  a  Mach  stem  at 
the  leading  edge  of  a  reflective  wall,  the  initial  value 
problem  for  the  wave  state  in  the  vicinity  of  the  caustic  of 
linear  theory  may  be  subject  to  the  development  of  a  wave 
jump  condition,  as  described  by  Peregrine  (1983) ,  and  the 
asymptotic  result  in  qualitative  similarity  to  the  Airy 
function  may  n.,t  develop. 

As  a  test  of  this  hypothesis,  two  wedge  geometries 
specified  by  6  =15  and  3 *25 3  were  tested  using  the 
parabolic  model  (6.2.9)  with  Pf  *0,  with  incident  wave 
steepnesses  of  k#Ao*0.1  and  0.2.  Waves  are  assumed  to  enter 
the  grid  parallel  to  the  x-axis,  and  the  y  boundaries  were 
assumed  to  be  totally  reflective,  with  one  boundary  located 
at  the  centerline  of  the  wedge-shaped  depression  and  the 
other  located  far  enough  away  from  the  wedge  that  waves 
reflected  at  the  caustic  do  not  reflect  back  into  the 
computational  domain  from  the  side  boundary  within  the 
x-distance  used.  Results  for  &  *15  and  6*25  are 
presented  in  comparison  to  results  obtained  by  the 
linearized  form  of  (6.2.9)  in  Figures  6.25  and  6.26, 


respectively,  where  normalized  amplitudes  perpendicular  to 
the  line  of  symmetry  y*0  are  plotted  for  k  x*40,  80,  and  120 


k  x  s  40 

O 


a)  k0A0  *  0.1 

Figure  6.26  Submerged  depression,  9  »  25°.  Normalized 
amplitude  vs.  distance  from  centerline  of 
depression 

a)  k  A  -  0.1;  b)  k^  -  0.2 

- linear  theory,  P^  »  0 

-  nonlinear  theory,  Pi  »  0 


[I] 


for  each  wave  steepness  tested.  For  the  case  of  a  *15  and 
kdA0  *0.1  (Figure  6.25a),  the  evolution  of  the  nonlinear 
wave  in  the  vicinity  of  the  caustic  is  not  qualitatively 
dissimilar  from  that  of  the  linear  wave.  The  maximum 
amplitude  near  the  caustic  is  seen  to  be  decreased  by  the 
effect  of  nonlinearity,  but  the  generation  of  a  reflected 
wave  and  a  large  amplitude  decay  in  the  shadow  zone  are 
apparent.  For  c*  15°  and  kdA0  *0.2,  however,  the  reflected 
wave  is  not  as  apparent,  and  a  broad  region  of  waves  which 
are  slightly  higher  than  the  incident  wave  develops  in  the 
vicinity  of  the  caustic.  Wave  amplitude  decays  much  more 
slowly  towards  the  center  of  the  wedge  in  the  area  of  the 
linear  shadow  zone. 


The  result  of  nonlinearity  is  even  more  accentuated 
in  the  results  for  Q  *25‘:,,  where  little  reflection  of  the 
incident  wave  is  apparent  for  either  incident  wave 
steepness.  In  order  to  accentuate  the  qualitative 
differences  between  the  linear  and  nonlinear  results,  plots 
of  the  instantaneous  surface  for  the  linear  wave  field  and 
for  the  nonlinear  wave  with  kftAfl  *0.2  are  given  in  Figures 
6.27  and  6.28,  respectively,  for  a  wedge  angle  of  6=25’. 


It  is  apparent  in  Figure  6.28  that  a  broad  wave  crest 
travelling  parallel  to  the  caustic  region  has  developed. 
Peregrine(1983)  argues  that,  in  the  event  that  the  jump 
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conditions  develop,  this  wave  crest  mast  continue  to  grow  in 
width  since  the  wave  energy  flux  continually  incident  on  the 
jump  boundary  cannot  be  balanced  by  the  flux  in  a  jump  of 
constant  width  and  height  unless  total  reflection  of  the 
incident  wave  occurs.  It  is  therefore  unlikely  that  the 
wave  field  will  evolve  into  its  asymptotic  state. 

Since  the  principle  conclusions  on  the  evolution  of 
the  nonlinear  wave  field  concern  the  absence  of  the 
development  of  a  reflected  wave,  the  results  for  the  wave 
steepness  k^A^  *0.2  and  both  d  *15°  and  &*25J  were 
recalculated  using  the  more  accurate  model  (6.2.9)  with 
Pj  *1/4  to  insure  that  reflections  were  not  being 
artificially  suppressed.  Although  some  increase  in  the 
magnitude  of  the  reflected  wave  were  noted,  the  results  for 
the  region  of  the  caustic,  where  nonlinear  effects  dominate. 


were  unaltered. 
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6.8.  Reflection  from  a  Caustic:  Current  Case 

As  a  first  example  of  waves  in  the  presence  of  a 
horizontally  varying  current/  we  consider  the  case  of  a 
linear  caustic  caused  by  reflection  by  a  following  current. 
The  results  are  then  qualitatively  similar  to  the  results  of 
section  6.7. 


A  stationary  current  distribution  is  established  and 
is  given  by 


U£«0 


O  j  Je.«  ^ 

o 


v  •  o  (i.n) 

Here,  tfU/jy  is  taken  to  be  a  constant  for  each  experimental 
case.  The  computational  domain  is  given  by 

0  £  ^  120 

Waves  are  incident  in  the  region  k#y<80  at  an  angle  v  to 

the  x  axis,  and  open  boundary  conditions  are  specified  on  — 

the  lateral  boundaries.  The  initial  wave  field  is  given  by 


*■*-  sL  - 


.-•.11  vv:  m. 


209 


A  {<>,•&)  =  J  A„c  j  b.(j4  » 

L  6  j  k?  » *» 

No  waves  are  present  initially  on  the  current.  The  physical 
situation  is  described  in  Figure  6.29,  and  represents  a 
straight  wave  maker  oriented  at  an  angle  }  to  the  y 
coordinate,  with  a  linearly  varying  current  distribution 
starting  at  the  end  of  the  wave  maker  and  increasing  in 
magnitude  with  distance  from  the  wave  maker.  The  depth  is 
taken  to  be  a  constant,  with  kfth=1.0. 


Snell's  law  applied  to  this  problem  requires  that 
the  x-component  X  of  the  wavenumber  vector  be  constant 
throughout  the  domain,  and  is  given  by 


J.  -  kt  z;-  &  =  cowifikt 


U.ts) 


The  wavenumber  is  then  uniquely  determined  throughout  the 
domain  according  to 


u>  =  a-  *  J.U. 


(t.  i  l) 


The  position  of  the  caustic  in  the  current  distribution  is 

f 

also  simply  determined,  and  occurs  where  k(y)*A.  Since  k 
then  equals  kflcos  i?  at  the  caustic,  the  dispersion  relation 


may  be  written  as 


CJ  =  &  huh  (  k0k  cos  *)\ 


'/i 


k^LLcAd 


or,  solving  for  (J, 


‘/t 


U.  =  &Q  -  f  )?a  7v»0^  ^£g£C>_)_| 


fc*  COi  & 


The  value  of  U  at  the  caustic  is  thus  uniquely  determined  by 
the  choice  of  Q  ,  which  gives  k^  in  turn. 


Two  sets  of  computations  were  run  with  £  *15 ’  and 
T*2.30sec.  for  each  case.  The  corresponding  value  of  k6  is 
lm“*  ,  and  U  at  the  caustic  is  calculated  from  (6.8.8)  to  be 

G(,  *  0.07 FT  n  /ttc 


Calculations  were  performed  using  the  numerical  scheme 
(6.2.11-12),  and  results  were  obtained  for  linear  waves  and 
for  values  k„A  *0.1  and  0.2,  as  in  section  6.7.  For  the 
first  set  of  calculations,  jU/jy  was  set  equal  to  0.05,  and 
the  caustic  occurs  at  k0y*81.51.  Results  are  presented  in 
Figure  6.30.  For  kftAo*0.1  (Figure  6.30a),  nonlinear  effects 
are  relatively  weak,  and  there  is  little  qualitative 


difference  between  linear  and  nonlinear  solutions. 


For 


Figure  6.30  Reflection  of  waves  from  caustic 
caused  by  a  following  current; 
3U/3y  =  0.05,  9  =  15° 

- -  nonlinear  solution; 

-  linear  solution 

a)  k  A  -  0.1;  b)  k  A  =  0.2 
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koAo*0.2  (Figure  6.30b),  a  broadened  wave  crest  develops  in 
approximately  the  same  position  as  the  first  antinode  in  the 
linear  wave  solution,  with  very  little  spreading  of  wave 
energy  beyond  the  caustic  being  apparent.  The  position  of 
the  caustic  is  indicated  by  the  vertical  line  in  the  figure. 

A  second  set  of  calculations  were  performed  with 
oU/«?y*0.025,  representing  a  more  slowly  varying  current 
distribution.  The  caustic  occurs  at  koy=83.02.  Results  are 
shown  in  Figure  6.31.  The  increased  effect  of  nonlinearity 
for  the  case  koAo»0.2  (Figure  6.31b)  is  dramatic.  A  broad 
wave  crest  develops  in  the  region  of  the  caustic  and  spreads 
significantly  beyond  the  caustic  position.  The  reflection 
of  waves  from  the  region  of  the  caustic  has  also  been 
suppressed  to  a  large  extent.  It  is  evident  that  a  wave 
jump  has  occurred,  and  that  the  conjugate  wave  state  is  able 
to  coexist  with  the  current  in  the  shadow  region  of  the 
caustic.  The  wave  would  be  expected  to  decay  with 
increasing  k0y  due  to  the  continual  increase  in  current 
velocity. 

A  possible  explanation  for  the  increase  in  nonlinear 
effects  with  decreasing  cU/->y  is  that,  for  progressively 
lower  current  gradients,  waves  entering  the  current  are 
refracted  more  slowly  and  thus  travel  a  greater  distance 


before  leaving  the  current  as  reflected  waves.  This 
increased  distance,  over  which  incident  waves  and  reflected 
waves  are  nearly  colinear  in  the  vicinity  of  the  caustic, 
allows  for  a  greater  span  of  nonlinear  interaction  between 
the  individual  waves  and  enhances  the  possibility  for  wave 
jump  conditions  to  develop.  An  analogous  behavior  would  be 
seen  in  the  results  of  section  6.7  if  the  slope  of  the  side 
walls  of  the  wedge  were  decreased. 


6.9 


Interaction  with  Rip  Currents 


Several  coastal  features  lead  to  the  presence  of  a 
narrow,  offshore-directed  jetlike  flow.  On  a  large  scale, 
ebb-tidal  flows  from  inlets  and  river  mouths  may  influence 
wave  propagation  over  large  stretches  of  adjacent  coastline. 
On  a  smaller  scale,  narrow,  intense  seaward  flows  known  as 
rip  currents  may  occur  on  otherwise  uniform  coastlines,  and 
lead  to  the  presence  of  on-offshore  channels  in  the  beach 
face  which  may  stabilize  the  rip  current  and  contribute  to 
considerable  offshore  sand  transport.  In  this  section,  we 
study  the  interaction  of  an  incident  wave  train  with  such  a 
current.  As  in  the  previous  section,  the  flow  pattern  of 
the  current  is  assumed  to  be  fixed. 

In  one  of  the  first  investigations  of  the  influence 
of  currents  on  results  of  the  ray  approximation  of 
refraction,  Arthur (1950)  studied  the  interaction  between  a 
normally  incident  wavetrain  on  a  plane  beach  with  a  flow 
field  designed  to  mimic  the  features  of  a  typical  rip 
current.  With  the  coordinate  x  oriented  offshore  and  the 
rip  centerline  located  along  x  and  at  y=0,  Arthur's  velocity 
field  is  given  by 

-(X/76.2)/4  -(£j/7.4Z)1/l 

U  3  o.oxi^s-  x  e  e 


-fr/7<.2)*/2 

V-  -  o.im  (z-  (x/74.2)  )  e 
•  er 

where  U  is  positive  in  the  offshore  direction  and  erf 
denotes  the  error  function.  Current  velocities  are  in 
meters/second.  The  bottom  topography  consists  of  a  plane 
beach, 

h  60  *  0.0 1  x  ,  (6.9.  z) 

Arthur  calculated  the  ray  pattern  for  an  incident  wave 
period  of  8  seconds;  his  results  are  shown  in  Figure  6.32 
for  comparison  with  the  refraction-diffraction  results 
discussed  below. 

Arthur's  solution  technique  involved  some  confusion 
over  the  relation  between  wave  rays  and  wave  orthogonals,  as 
evidenced  by  the  straight  "rays"  to  either  side  of  the 
narrow  rip;  the  ray  paths  should  be  convected  towards  the 
centerline  of  the  rip  by  the  longshore  current  perpendicular 
to  the  wave  orthogonals.  Christof fersen(1982)  noted  this 
discrepancy  during  his  development  of  a  ray-tracing 
procedure  for  waves  in  the  presence  of  a  current,  and 
Jonsson,  Chr istof f ersen  and  Skovgaard  (1983) ,  using 
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Chr istoffersen* s  procedure,  have  presented  a  corrected 
version  of  Arthur's  results  which  also  include  estimates  of 
the  local  wave  amplification  with  respect  to  the  incident 
amplitude  at  h»100m.  Jonsson,  Christof fersen  and 
Skovgaard's  results  suggest  a  decrease  in  wave  height  in  the 
center  of  the  rip;  this  appears  to  be  a  result  of  the 
termination  of  several  rays  and  a  subsequent  underestimation 
of  the  local  wave  energy  density.  These  results  point  out 
the  difficulties  of  interpreting  refraction  patterns  in 
regions  of  multiple  ray  crossings. 

In  this  section,  we  use  the  parabolic  model  given  by 
(6.2.11-12)  to  investigate  the  interaction  between  a 
normally  incident  plane  wave  and  a  fixed  current  pattern 
given  oy  (6.9.1).  The  deviation  of  the  free  surface  due  to 
the  0(1)  current  is  neglected. 

At  large  distances  from  the  centerline  of  the  rip 
current,  the  wavefield  is  characterised  by  normally  incident 
waves  and  a  longshore  current  distribution  with  no 
on-offshore  component.  Consequently,  the  wave  crests  and 
snoaling  characterisitcs  should  not  be  affected  by  the 
current  at  lateral  boundaries  far  from  the  rip.  In  order  to 
place  the  lateral  boundaries  close  to  the  rip  and  reduce  the 
size  of  the  computational  domain,  the  lateral  boundary 
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conditions  are  replaced  by  a  one-dimensional  model  for 
shoaling  waves,  given  by 


A* -j.(k-K')*  *JL  C.XA 

x  2C^ 


•f  Jj/  /4  +  J?,  DM  I  A  * 


j) 


i;: 


which  is  just  the  one-dimensional  counterpart  of  (4.4.4). 

Breaking  is  included  in  all  calculations  in  this  section.  ’-j 

We  first  consider  linear  theory.  As  mentioned  in  18 

section  4.4,  the  parabolic  approximation  for  wave-current  •*-; 

interaction  derived  by  the  splitting  method  contains  errors 
in  y-derivative  terms.  We  first  compare  results  obtained  ISO 

using  the  approximate  equation  (4.4.33),  written  in 
finite-difference  form  as  (6.2.11),  to  results  obtained 
using  the  correct  equation  (4.4.7),  written  in  a  similar  *£ 

finite-difference  form.  The  computational  grid  is  given  by 

Osr  x'  &  ’T7.  a"  m  «  ~  87. 1*  s  s  ^ 


with  x'*0  located  offshore  and  the  shoreline  located  at 
x'«300m.  The  values  of  ax  and  Ay  are  2.5m  and  1.5m, 
respectively.  A-  plot  of  the  instantaneous  surface  elevation 
for  an  amplitude  A0(h*6m)  *0.1m  and  a  wave  period  T=8sec.  is 
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given  in  Figure  6.33.  Note  that  the  wave  amplitude  is 
strongly  amplified  along  the  centerline  of  the  rip,  and  that 
discontinuities  develop  in  the  wave  phase  surface,  leading 
to  the  presence  of  nodes  in  wave  amplitude  along  the  lateral 
boundaries  of  the  core  of  the  offshore  directed  rip.  The 
amplification  seen  here  is  in  direct  contradiction  to  the 
results  presented  by  Jonsson,  Christof fersen  and 
Skovgaard (1983) .  Wave  breaking  occurs  on  the  rip  centerline 
at  x'»240m,  or  60m  from  the  shoreline. 

The  region  of  strong  amplification  and  large 
currents  along  the  centerline  of  the  rip  is  likely  to  be  a 
good  indicator  of  the  amount  of  error  involved  in  using  the 
approximate  equation  (4.4.33).  Linear  wave  results  were 
obtained  for  an  initial  amplitude  of  A#=0.01m  and  T=8sec. 
using  both  the  exact  (4.4.7)  and  approximate  (4.4.33) 
equations.  For  this  amplitude,  breaking  is  confined  to  the 
last  shoreward  grid  row.  Plots  of  amplitude  relative  to  the 
incident  wave  are  given  in  Figure  6.34  for  transects  y=0 
along  the  rip  centerline  and  for  values  of  x' *197. 5m, 
222.5m,  247.5m  and  272.5m.  The  results  show  a  decrease  in 
maximum  amplitude  in  the  approximate  solution  relative  to 
the  results  of  the  exact  equation  in  the  region  of  strong 
focussing  near  tne  shoreward  end  of  the  rip,  indicating  a 
loss  in  wave  action  density  due  to  tne  errors  in  terms 


Figure  6.34  Amplitude  relative  to  incident  wave  for  waves 
interacting  with  rip  current.  Comparison  of 
results  of  exact  equation  (4.4.4)  and  approxi¬ 
mate  equation  (4.4.33)  :  linear  waves  - 

(4.4.4) ;  -  (4.4.33) . 
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involving  the  y  velocity  component  V  in  (4.4.33).  These 
results  indicate  that  it  is  not  desirable  to  use  the 
approximate  form  (4.4.33)  if  detailed  estimates  of  local 
wave  amplitude  are  required.  The  transects  for  constant  x' 
values  clearly  show  the  development  of  the  node  in  the 
amplitude  as  the  discontinuity  in  wave  crests  develops. 

Due  to  the  strong  wave  focussing  in  the  area  of  the 
rip  current,  nonlinear  effects  are  also  likely  to  be 
important  in  this  example.  As  in  several  of  the  cases 
studied  in  section  6.6,  nonlinear  effects  occur  in  this 
example  at  ursell  numbers  too  large  for  Stokes  wave  theory 
to  be  valid,  if  reasonable  values  of  initial  wave  amplitude 
are  chosen.  In  order  to  investigate  the  effect  of 
nonlinearity  in  this  example,  we  incorporate  the  effect  of 
amplitude  dispersion  in  the  model  using  the  empirical 
dispersion  relation  of  Walker(1976) .  The  effect  of 
nonlinearity  is  included  in  the  model  by  using  the 
linearised  governing  equation  but  calculating  k  based  on  the 
revised  dispersion  relation,  as  suggested  in  section  3.5. 

The  correction  to  linear  wave  theory  is  given  by 
Walker  as  a  modification  to  the  linear  phase  speed. 
Walker's  phase  speed  Cv  is  given  by 


in  the  present  notation/  where  C  is  the  phase  speed 
determined  by  the  linear  dispersion  relation.  Walker 
obtained  this  one  parameter  fit  by  comparison  with  data 
obtained  from  laboratory  waves  shoaling  on  a  plane  beach.  A 
dispersion  relation  may  be  determined  from  (6.9.4),  and  is 
given  by 


W l-  <3  te„  *  i  lJr  )  ^  f ( 1  *  !4L) 


Walker's  modification  is  purely  empirical  and  it  is 
desirable  to  gain  some  insight  on  how  his  results  compare  to 
theory,  in  order  to  provide  this  information,  we  compare 
results  for  Walker's  phase  speed  prediction  to  values  given 
by  the  stream  function  wave  theory  of  Dean(1965) ,  which  is 
valid  for  all  water  depths.  Values  are  compared  for  a  range 
of  local  wave  heights  relative  to  the  local  breaking  wave 
height  as  determined  by  the  stream  function  theory,  and  for 
a  range  of  relative  depths  h/Lft ,  where  L0  is  the  linear  deep 
water  wavelength  given  by 
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Values  of  H(*2|Aj)/Hb  Of  0.25,  0.5,  0.75  and  1.0  were  used. 
Stream  function  values  were  obtained  from  the  tables  of 
Dean(1974) .  Results  of  Walker's  modified  phase  speed  are 
presented  as  a  %  error  relative  to  stream  function  theory  in 
Figure  6.35.  Linear  theory  results  are  included  for 
comparison.  The  range  of  h/L0  values  relevant  to  the 
present  problem  extends  from  0.060  for  the  offshore  boundary 
to  0.012  near  x'=240,  where  breaking  occurs  in  the  linear 
theory  for  the  choice  Ao»0.1ra.  it  is  apparent  from  Figure 
6.35  that  Walker's  modification  provides  a  significant 
improvement  over  linear  theory  except  for  the  region  close 
to  the  offshore  boundary,  where  wave  amplitude  is  small. 
However,  for  small  amplitudes  the  relative  error  in  either 
Walker's  or  linear  results  are  on  the  order'  of  2%.  It 
appears  that  Walker's  empirical  formula  provides  an  adequate 
estimate  of  the  phase  speed  even  for  steep  waves  near 
breaking.  We  remark  that  it  would  be  desirable  to  obtain 
more  detailed  forms  of  the  empirical  dispersion  relation 
designed  to  cover  the  entire  range  of  depth  values,  since 
Walker's  relation  is  suitable  mainly  in  shallow  water  and 
approaches  linear  theory  asymptotically  in  deep  water. 

For  the  present  application  to  wave-current 

interactioa,  we  modify  Walker's  dispersion  relation 


according  to 


Error 


.002  .005  .01  .02 


.05  .1  .2 


.5  1.0 


« 

Error 


.002  .005  .01  .02  .05  .1 


.5  1.0 


Figure  6.35  Error  in  predicted  phase  speed  relative 
to  results  of  stream  function  theory, 
a)  H/Hb  *  0.25,  b)  H/Hb  *  0.5, 

c)  H/Hb  *  0.75,  d)  H/Hb  *  1.0 

-  Walker  (1976) ;  -  linear  theory 
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CJ  -  C*  ^  t*/  ^  =  *  few  U- 


i.1.7 


wnere 


The  wavefield  for  the  case  A  =0.1m  and  T=8sec.  is  shown  in 

O 

Figure  6.36.  Plots  of  amplitude  for  the  same  transects  as 

in  Figure  6.34  are  given  in  Figure  6.37  in  comparison  to 

« 

linear  results.  The  transect  at  x' *272. 5m  is  inside  the 
surfzone  for  both  the  linear  and  nonlinear  results. 


It  is  apparent  from  Figure  6.37  as  well  as  from  a 
comparison  of  Figures  6.36  and  6.33  that  the  effect  of 
nonlinearity  due  to  the  empirical  shallow  water  dispersion 
relation  is  relatively  subdued  in  comparison  to  the  striking 
examples  of  sections  6.7  and  6.8.  An  invalid  application  of 
Stokes  theory  to  the  present  example  leads  to  the  prediction 
of  strong  wave  jump  conditions  in  the  region  of  the  rip; 
this  effect  is  absent  from  the  present  calculations.  The 
nonlinear  effects  are  limited  to  a  reduction  in  wave 
amplitude  along  the  centerline  of  the  rip  and  a  slight 
broadening  of  the  region  of  focussed  waves  as  evidenced  by 
the  position  of  the  partial  nodes  in  the  amplitude 
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Amplitude  relative  to  incident  wave  for 
waves  interacting  with  rip  current. 
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transects.  The  results  are  similar  to  the  nonlinear 
modification  seen  in  the  data  for  waves  passing  over  the 
elliptic  shoal  in  section  6.5. 

The  results  of  this  section  indicate  that  it  may  be 
possible  to  make  use  of  a  variety  of  dispersion  relations 
within  the  context  of  the  linear  wave  model  in  order  to 
model  phenomena  beyond  the  limits  of  validity  of  Stokes 
theory.  Before  any  conclusions  regarding  this  possibility 
can  be  drawn,  it  will  be  necessary  to  investigate  the 
behavior  of  the  models  in  comparison  to  existing  or  future 
data  for  waves  propagating  and  shoaling  in  shallow  water. 

The  benefits  to  be  achieved  in  the  event  of  the  validity  of 
thisr  approacn  are  clear,  since  existing  numerical  approaches 
depend  on  the  use  of  Boussinesq  or  Kor teweg-deVr ies 
equations,  which  require  more  involved  computational 
procedures  than  the  parabolic  equation  method  used  here.  P! 


Ut 


CHAPTER  7.  SUMMARY  AND  SUGGESTIONS  FOR 


FURTHER  RESEARCH 


In  this  study,  we  have  sought  to  provide  a  model  for 
progressive  surface  waves  in  the  presence  of  a  large, 
horizontally  varying  current.  The  formulation  includes  the 
lowest  order  effects  of  nonlinearity  in  a  manner  consistant 
with  the  Stokes (1847)  expansion  in  a  small  parameter  based 
on  the  wave  steepness.  In  addition,  the  depth  is  allowed  to 
vary  somewhat  more  rapidly  than  in  the  traditional  mild 
slope  approximation,  and  terms  proportional  to  the  square 
and  derivative  of  the  local  bottom  slope  are  retained. 

In  Chapter  2,  we  obtain  a  solution  for  the  velocity 
potential  for  a  plane  progressive  wave  to  OCe1)  in  wave 
steepness  and  0(2)  in  bottom  slope.  Then,  using  the 
variational  principle  of  Luke(1967),  we  obtain  a  Lagrangian 
governing  the  entire  fluid  motion  in  the  absence  of  viscous 
effects,  given  in  Chapter  3.  Variations  of  the  unaveraged 
form  of  the  Lagrangian  with  respect  to  the  unknown  dependent 
variables,  following  a  perturbation  scheme  developed  here. 
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yields  equations  governing  the  instantaneous  fluid  motion. 

The  lowest  order  problem  reproduces  the  nonlinear 
shallow  water  theory  governing  the  0(1)  mean  motion.  This 
result  could  be  further  expanded  for  the  case  of  small  bft /h^ 
to  yield  the  Boussinesq  equations. 

For  the  leading  order  wave  motion  (0(c))/  we  obtain 

,  O'  r 

an  equation  including  terms  to  O(i')  in  and  ,  governing 

/V 

the  potential  .  Invoking  the  simplifying  assumption  of 

plane  wave  motion  reduces  the  nonlinear  contribution  to  a 

coefficient  proportional  to  the  square  of  the  local 

O' 

amplitude  multiplying  the  potential  <^>  .  This  term  has  the 
effect  of  modelling  amplitude  dispersion  arising  at  0 (£  ). 
Interaction  with  the  mean  0(  2 *)  wave-induced  motion  is  also 
included . 

x 

At  0(  s  ) ,  a  set  of  forced  equations  representing  a 
free  surface  boundary  condition  and  continuity  equation  for 
the  entire  0( O  motion  is  obtained.  Taking  the  mean  over  a 
wave  period  yields  the  continuity  equation  for  wave  wave 
induced  mean  flow,  as  found  by  Whitham (1962) ,  and  a  forced 
wave  equation  governing  the  0(cl)  wave-induced  flow  is 
obtained.  Retaining  the  the  full  equations  at  0(i‘)  yields 
equations  which  could  be  used  to  calculate  the  forced  second 


harmonic  of  the  fundamental  wave  g?  • 


The  resulting  formulation  is  a  correct 
representation  of  the  governing  equations  for  plane 
progressive  waves  in  the  absence  of  strong  reflection.  This 
is  demonstrated  in  Chapter  4,  where  we  show  the 
correspondence  between  the  second  order  wave  equation 
formulation  and  the  nonlinear  Schrodinger  equation  governing 
the  evolution  of  the  amplitude  envelope.  This 
correspondence  points  out  the  potential  capability  for 
modelling  wavefields  with  narrow-banded  frequency  spectra 
and  an  identifiable  carrier  frequency.  However,  in  the 
remainder  of  this  study,  we  choose  to  limit  our  attention  to 
steady,  monochromatic  wavefields.  parabolic  approximations 
for  waves  in  the  mild  slope  formulation  are  derived  in 
Chapter  4  and  used  to  study  a  number  of  examples  in  Chapter 
6.  A  detailed  comparison  with  laboratory  data  for  the  case 
of  waves  propagating  over  a  submerged  shoal  demonstrates  the 
validity  of  the  parabolic  equation  method  and  indicates  that 
accurate  solutions  may  be  obtained  using  this 
computationally  efficient  method.  This  result  has 
significant  bearing  on  the  potential  for  development  of  wave 
models  for  engineering  application. 

In  contrast  to  the  case  for  progressive  waves,  the 
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formulation  contains  basic  ambiguities  in  the  case  of 
arbitrary  wave  motions,  both  in  terms  of  the  evaluation  of 
several  of  the  second  order  terms  in  the  bottom  slope  as 
well  as  in  the  evaluation  of  the  wavenumber  for  the  case  of 
waves  on  a  current.  Both  evaluations  require  that  the 
direction  of  a  wavenumber  vector  be  specified  relative  to 
the  direction  of  the  current  vector  or  vector  gradients  of 
various  quantities  in  the  physical  domain.  The  model  does 
not  directly  calculate  the  local  wave  direction,  and,  in  the 
case  of  partial  standing  waves,  this  quantity  is  not  well 
defined.  The  problems  associated  with  this  ambiguity  are 
pointed  out  in  Chapter  5,  where  we  study  the  reflection  of 
waves  by  steep  depth  transitions. 

One  path  along  which  future  research  could  proceed 
is  towards  removing  some  of  the  limitations  of  the  present 

I 

theory.  An  expansion  of  y  without  assuming  a  plane  wave 
form  of  the  solution  a  priori  may  alleviate  the  difficulties 
encountered  in  Chapter  5.  Further,  the  limitations  imposed 
by  the  assumed  irrotationality  of  the  motion  should  be 
investigated  and,  if  possible,  eliminated,  since  it  is  clear 
that  most  physical  cases  of  interest  will  involve  currents 
with  rotational  velocity  distributions.  Finally,  although 
the  effect  of  viscosity  is  neglected,  we  have  shown  that  it 
is  possible  to  effectively  parameterize  viscous  effects  in 
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the  model  simply  through  imposing  the  condition  of  a  complex 
frequency  or  wavenumber.  The  resulting  equations  are  able 
to  effectively  model  both  gradual  and  strong  dissipative 
effects  including  wave  breaking.  For  applications  to  large 
scale  problems,  it  would  also  be  desirable  to  include  the 
facility  to  allow  for  wind  generation  of  waves  as  well  as 
deep  water  wave  breaking  due  to  excessive  wave  steepness. 

The  parabolic  and  Sc^rodinger  equation  wave  models 
as  formulated  here  are  suitable  for  direct  use  in  models  of 
wave  induced  circulation  in  the  nearshore  region.  Most 
currently  existing  models  such  as  those  by  Wu(1983)  and 
Ebersole  and  Dalrymple  (1980,  see  also  Kirby  and 
Dalrymple (1982)  for  a  recent  summary)  use  a 
finite-difference  refraction  scheme  and  solution  of  the  wave 
energy  equation  to  provide  the  local  wave  amplitude  and  wave 
angle  for  use  in  the  momentum  equations  for  the  mean  flow. 
These  refraction  schemes  do  not  allow  directly  for 
diffraction  and  do  not  provide  detailed  information  in 
regions  of  strong  focussing  and  ray  crossing.  The  models 
should  therefore  be  especially  useful  in  applications  to 
regions  with  strong  currents  such  as  rip  currents  and  ebb 
tidal  flows  from  inlets. 

The  limitation  of  the  validity  of  Stoke' s  theory  in 
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shallow  water  places  a  severe  constraint  on  the  general 
utility  of  the  models  formulated  here.  Most  real  waves  in 
the  nearshore  zone  are  more  appropriately  modelled  by  the 
Boussinesq  equations.  However,  it  has  been  shown  in  section 
6.9  that  models  of  a  more  empirical  nature  may  be  developed 
by  the  direct  substitution  of  any  given  dispersion  relation 
into  the  linear  form  of  the  governing  equations.  It  may  be 
possible  to  develop  methods  by  which  the  basic  linear 
propagation  model  may  be  made  to  mimic  any  desired  real  wave 
process  (such  as  the  propagation  of  broken  waves  in  the 
surfzone)  by  a  suitable  choice  of  dispersion  relation  and 
dissipative  mechanism.  It  is  not  clear,  however,  that  such 
a  model  would  be  capable  of  maintaining  a  proper  balance 
between  relevant  quantities  such  as  amplitude,  energy  flux 
and  phase  speed. 
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APPENDIX  A.  THE  STOKES  SOLUTION  TO  0 ( £  )  IN  THE 


PRESENCE  OF  A  CURRENT 


In  this  appendix,  the  series  solution  of  Chu  and 
Mei (1970a)  is  extended  to  include  the  case  of  an  0(1)  mean 
current  JJ(x,t)  ,  given  by 

where  <|>-|o(x,z,t)  is  the  leading  term  in  a  Stokes  expansion 
and  is  0(  6  1  )  .  The  term  is  assumed  to  be  slowly  varying 
in  time.  The  ambient  current  is  then  of  0(6  '£  ),  or, 
equivalently,  0(1)  in  the  present  context.  In  this 
appendix,  the  notation  of  Chu  and  Mei (1970a)  is  used  for 
convenience,  with  the  consequence  that  the  ordering 
subscripts  for  quasi-steady  terms  and  0(  £  )  contributions 
are  respectively  one  order  smaller  or  larger  than  in  the 
main  body  of  the  text. 

The  series  solution  of  Chu  and  Mei  for  (j>  and  is 
given  by  (after  substituting  the  modulation  scale  S  for  the 


equivalently  valued  parameter  £  where  appropriate) 
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and  *  denotes  the  complex  conjugate  of  the  preceeding  term. 
The  terms  fJo  -f  ZQ  are  given  in  (2.3.2)  and  (2.3.5).  Here, 
the  '  s  are  given  in  terms  of  the  absolute  frequency  CO  : 


The  lowest  order  quasisteady  term  in  the  expression  for  (fi 

i 

is  ^0,  which  is  related  to  the  second  order  wave- induced 

current.  It  is  desired  to  alter  the  expansion  to  include  a 

-f 

term  of  0(  c  )  which  leads  to  an  0(1)  current  U  upon 


differentiation. 


For  the  standard,  slowly  modulated  form  of  the 
Stokes  expansion,  neglecting  <^(  and  ,  it  is  well  known 
that  the  inclusion  of  a  constant  mean  current  U,  leads  to  an 
expression  for  the  relative  (or  intrinsic  )  frequency  CT  , 

$-  CO  -  64-0 


The  terms  in  the  expansion  of  <p  are  merely  adjusted  to 
account  for  the  frequency  shift  by  means  of  a  Galelean 
transformation;  see  Thomas  (1979 )  .  Thus  <£>(|  and  are 

altered  to  the  form 
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r 


n  Jz 


20 


We  must  now  determine  the  corresponding  form  of  the  fast 
modulation  term  ,  as  well  as  any  higher  order  terms  due 
to  the  addition  of  the  slowly  varying  term  at  0(  6  1  ) .  As 
mentioned  in  Chu  and  Mei (1970a) ,  we  may  obtain  this 
expression  in  the  context  of  an  expansion  for  linear  waves; 
this  is  particularly  clear  if  we  consider  as  a  term  of 
0(€t,'6).  Correspondingly,  we  choose  expansions  of  the  form 
(  see  Chu  and  Mei (1970a)) 

€*<E  <l  z  (a.u) 


where  jm  {  <1  since  we  do  not  need  to  consider  the  higher 

harmonics  arising  in  the  nonlinear  problem.  The  governing 
equations  are  written  in  the  form 

<Pxx  -  =  °  CM*) 

3  (A.< 

after  elimination  of  from  the  free  surface  boundary 
conditions.  Equations  (A.9a-b)  are  sufficient  to  obtain  the 
results  desired  here.  The  free  surface  boundary  condition 
is  used  to  provide  consistancy  conditions  which  will  be 
obtained  separately  by  the  variational  approach  in  Chapter 
3.  As  an  exception,  we  will  use  (A. 9c)  to  verify  the 
dispersion  relation  (2.1.6). 

The  expansions  (A. 8)  are  substituted  into  the 
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equations  (A. 9),  leading  to  the  set  of  equations 

j  (A.iod 

‘  F'~  i  **"A*  (AM) 

=  <s'**-  j  a*1”  Ca.idc.) 

where,  for  each  order  n,  the  terms  R,  G,  and  F  are 
determined  in  terms  of  lower  order  quantities  or  corrections 
to  the  absolute  frequency  Uy  .  The  Taylor  series  expansion 
of  the  free  surface  condition  (A.9.c)  is  made  with  respect 
to  the  0(1)  free  surface  position  z=  f^00,  since  this  is  the 
total  depth  occupied  by  the  current  <jL,0  in  the  absence  of 
waves.  We  define  a  shifted  vertical  coordinate  z’*z -  J 
the  free  surface  conditions  are  then  applied  to  the  position 
z '  =0 ,  while  the  total  depth  influencing  the  wave  motions  is 
given  by  h=ho  +  .  We  find  that,  to  lowest  order, 

R  =  f  =  a  =  o 

and  (A.IOa-b)  are  satisfied  by  ^_l6(x,y,t)  .  The  lowest 

order  componant  of  the  ambient  current  thus  does  not  vary 


over  the  depth,  which  is  consistent  with  the  neglect  of 
friction  in  the  formulation. 


Turning  now  to  the  problem  for  m=0,  n=0,  it  is  found 

that 

\ 

leading  to  a  homogeneous  problem  for  i^.  Foda  and  Mei(1981) 

have  used  <k  to  allow  for  the  resonant  growth  of  a  long 
T  00 

forced  wave  driven  by  the  0(  c,  )  wave  field.  This  effect  is 

I 

not  central  to  our  present  concerns,  and  we  neglect  <fao  . 
The  problem  at  m=Q,  n=l  is  altered  from  Chu  and  Mei's 
problem  due  to  forcing  by  the  n=«-l  current,  and  the 
coefficients  R,  ,  F.  ,  and  GIa  take  on  the  values 

lo  10 

R„* 

H,  =  K  ■  v*  <k-u 

and 

whereas  in  Chu  and  Mei's  study  the  corresponding  problem  is 

l 

homogeneous.  The  inhomogeneous  problem  for  {>l0  is  given  by 


-b  <  z!  cA.i/a) 


and 
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^lO  -  ^6  (A.  (It) 

In  previous  studies  where  the  largest  current  is 
smaller  than  0{1) ,  has  been  identifiable  completely  as 
the  wave-induced  current  component.  Here,  the  homogeneous 
solution  of  (A. 13)  ,  denoted  as  ,  will  be  identified  with 

the  wave  induced  current.  The  remaining  problem  for  the 
particular  solution  can  then  be  solved  directly  using 

(A. 11  d,c),  leading  to  the  result 

(j>"  a  -  Ck+~l')  CA.Il) 

The  potential  <j?h  represents  a  higher  order 
correction  to  the  0(1)  ambient  current,  and  could  play  a 
similar  role  in  the  overall  wave-current  problem  as  the 
wave-induced  current  ,  in  the  sense  that  both  may  lead  to 

small  wavenumber  corrections  at  third  order  through  wave 

\U 

current  interaction.  It  should  be  noted  that  q>l0  is 
essentially  an  0(  S  )  correction  to  ,  as  is  ,  and  in 
that  sense  can  also  be  dropped  from  consideration  in  the 
context  of  a  mild-siope  formulation. 

Turning  to  the  problem  for  the  wave  component  at 
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n*l,  m=l,  we  find  that  R|(  *Fi(  *0,  while 

6„  =  {iu lz-vk^'  * 

Using  (A. 6)  as  a  definition  for  CT  ,  the  free  surface 
boundary  condition  (A. 10c)  can  be  written  as 

-  ffL4u  5=0  » 

The  solution  to  the  homogeneous  problem  (A.IOa-b) 
by 

4u  *  "44* 

La¬ 
in  agreement  with  (A. 7) ,  and  the  value  of  the  coefficient 
has  been  chosen  to  give  an  amplitude  A([  for  the  linear  wave  *v 

i 

component.  Substituting  £p)(  into  the  free  surface  boundary  m 

condition  (A. 13)  gives  the  result 

,Sj 

thus  verifying  the  dispersion  relation  given  in  section 

(2.1).  d 

♦  •/. 

At  n*2,  m*l,  and  RtJ  take  on  the  values 

r-*j 

to* 


Ml 


a 

The  solution  for  is  then  given  by 


r  /•  ? 

<j>w  »  cos/j  c)  --j^jR-u  Stilt  k &+%)</.%  r  + 

^  -A 

+  T/b/f  k(h  +  i')4^  C-l  +  ~g£j  cos/i  k(h+^)d^  (AM] 


where  and  are  coefficients  of  the  homogeneous 
solution,  and  where  we  have  used  the  method  of  variation  of 
parameters  (see,  for  example,  Hildebrand (1976) ,  p.24-26). 
Application  of  the  bottom  boundary  condition  (A. 10b)  gives 


the  result 


ci  *  Fv  /  k 


while  Cj  remains  formally  arbitrary.  After  performing  the 
required  integration  and  choosing  a  value  for  c(  consistant 

»  I 

with  the  form  of  w  ,  ytJ  is  given  by 
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<K,  =  -A%Av±io  ~2^U  ^  °i,  k(k*i')  flt  + 

+  ujdk+i')  suMthad  + 

e,,k  kk  (A.  is) 


Neglecting  the  homogeneous  tern  in  A  leads  to  the  solution 
of  Biesel (1951) ,  corrected  for  the  presence  of  a  current. 
Chu  and  Mei  (1970a)  noted  that  <|ja  becomes  unbounded  in  the 
limit  of  deep  water.  Therefore  A,j  is  retained  and  chosen 

I 

so  that  yu  remains  formally  bounded  in  the  deep  water  limit, 
yielding  the  results  of  Chu  and  Mei  corrected  for  the 


3 


presence  of  a  current: 


< fill  ~  /"  *  *  ^rj^nj 


(A.  It.) 


where  is  given  by 


N  t 


i  a 


The  appropriate  form  of  the  wavelike  componants  of  Chu  and 
Mei's  solution  is  thus  obtained  for  our  study  simply  by 


=  1 


replacing  the  absolute  frequency  CO  by  the  relative 
frequency  C~  in  each  term  in  (jj>  .  It  should  be  noted  that 
will  in  general  be  a  function  of  the  variation  of  the 
0(1)  current,  thus  greatly  complicating  the  expressions 
involving  and  derivatives  of 


The  results  of  the  perturbation  expansion  give  the 
form  of  and  r?  required  for  the  problem  formulated  in 
Chapter  2.  They  are  given  by 


i  *-/  ,  .A  \  A 

<j>  -  S  <p_l0  +  e<fu  e  +  *  +  fLI  A  +  * 
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where  the  two  values  of  ^  correspond  to  the  split  of  • 
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The  Lagrangian  L  for  an  expanded  dependent  variable 

4)1 

£■  €  f*.  /  n»0,l,2,...,  may  be  derived  by  substituting  the 

assumed  expansion  for  f  into  the  appropriate  form  of  L.  By 
retaining  the  ordering  parameter  £  »  L  may  then  itself  be 

expanded  in  the  form 

L  a  j  O  (B.L) 

If  L  is  in  general  nonlinear ,  then  lowest  order  terms  which 
are  quadratic  in  f^  will  appear  first  in  L  ;  subsequently, 
variation  of  with  respect  to  f^  leads  to  a  linear 

governing  equation  for  fw,  containing  forcing  terms 
involving  lower  order  terms  fK,  m<n.  The  exception  to  this 
rule  may  occur  at  0(1) ,  or  n*0,  where  the  governing  equation 
for  an  9(1)  quantity  fQ  may  not  be  reduced  to  a  linear  form; 
this  will  be  seen  in  the  problem  for  the  ambient  current  in 
Chapter  3. 

As  an  example,  consider  the  one-dimensional  Duffing 
equation : 


f  +  f  +  ef  -  o  (g.i) 

A  Lagrangian  corresponding  to  (B.2)  may  be  written  as 

L  •  X  -  _£  +  eX*  (m) 


/.*  X  -  X  +  <=X  (?•■*) 

z  i  v 

Here,  (  )  denotes  differentiation  with  respect  to  the 

independent  variable.  Let  an  expansion  for  f  be  given  by 

f  *  6*^  5^20  .  (B.v) 

(B.4)  is  then  substituted  into  (B.3)  and  ordered  into  the 


form  (B.l) ,  where 


hi  ~  ££  +  il/t 


(h.Sk) 

(&.  ?y) 


X  +  XX  *  ££-  £  £  (s-^) 


Proceeding  in  a  consistant  manner,  the  lowest  order 
contribution  to  L,  given  by  L0  ,  is  varied  with  respect  to 
the  lowest  order  unknown  contribution  to  f,  given  by  f c  , 
yielding,  after  partial  integration 


Then,  regarding  as  known  (as  determined  by  (B.6a)  ),  the 
next  lowest  order  term  Lj  is  varied  with  respect  to  the 
lowest  order  unknown  f|  ,  yielding  (B.6a)  as  a  redundant 
condition.  Thus  L t  ,  which  is  linear  in  the  lowest  order 
unknown,  yields  no  information  concerning  the  unknown. 
Then,  L £  is  varied  by  the  still  unknown  f(  ,  yielding  the 
governing  equation 

tea) 

This  process  may  be  continued  indefinitely.  (B.6a)  and 
(B.6b)  are  identical  to  the  results  obtained  by  applying  a 
regular  perturbation  scheme  directly  to  (B.2) ,  and  (B.6b)  is 
known  to  produce  a  secular  result,  which  may  be  removed  by  a 
stretching  of  the  independent  variable.  For  the  water  wave 
problem  to  be  studied,  the  corresponding  secular ity  appears 
at  0(d*),  resulting  from  ;  since  it  will  be  sufficient 
to  carry  the  calculations  to  L  ^  ,  the  perturbation  scheme 
will  be  regular  in  nature. 

While  the  treatment  presented  here  for  the  Duffing 
equation  does  not  represent  a  computational  advantage  over 
the  normal  procedure  of  perturbing  the  governing 
differential  equation  directly,  the  advantage  becomes 
apparent  in  the  case  of  problems  containing  cross-spaces, 
where  the  form  of  f  is  known  and  separable.  Then,  use  of 


the  Lagrangian  integrated  over  the  cross  space  leads 
directly  to  variational  principles  governing  motion  in  the 
propagation  space.  The  corresponding  governing  equations 
are  obtained  in  the  multiple  scale  expansions  as  solvability 
conditions  related  to  higher  order  inhomogeneous  problems. 


»** « *-»?*-■ 


■>. . y~ 


ml 


APPENDIX  C.  DERIVATION  OF  L  AND  THE  EULER  EQUATION 


Due  to  the  complexity  of  the  calculations  leading  to 
the  expression  for  L  and  the  resulting  Euler  equations,  the 
derivations  are  presented  here  in  some  detail,  with  L  being 
derived  in  section  1.  In  section  2,  the  Euler  equation  for 
the  wave  component  is  derived. 


C.l  The  Lagrangian  L 


The  series  form  for  ^  and  are  given  by 

(j)  a  S"'<£0  ■+  e{jlt  -A.%(ux{u  +  fa  *  4>  + 

+  S  •+  6l  \)"t  +  0(&\  S  ) 

(C.  I.  I ) 

and 

+  fe«2«  +  o45)  (£. /.  0 

Also,  is  expanded  in  powers  of  S  to  give 


where 

4>1  *  6='W) 

and 

<h"  =  -Lktii  vfk'  -  Ch*-i)\ h -vt  4' 

*  (C.I.S-) 

The  evaluation  of  L  requires  the  quantities  <j>^  and  V<^  , 

which  are  given  by 

<j>t  *  $,[  +  S* *  6 

+  +  -  X>  - 

4  (C./.4) 

where  we  have  invoked  the  condition  that  d/^~  of  the  ambient 
current  is  OfS*),  and 

W  = 

1  ^ 

’  V»  &'  «•  £V,  4"  ♦  4  S  V,  }  f.  iZeijfii  1  k, 

+  6  {  ^  S  4,  j-'j  1  + 


in 
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+  +  6f  j-io lz  h  + 

+  6  4o  2  4  ]  4-  a 


(C.f.7) 


The  expression  for  (V^)1*  only  needs  to  retain  terms  to 

■.  rr  V  .,  / 

0(a*,  a  )  and  0(d),  and  can  then  be  written  as  (  where  <j  s 


have  been  dropped  for  convenience) 


(s£i  3  +  'Vhio-Vhf*  +  z'-fafS ^.z 

z  z  7  X 

+  fejLf 

2.  Z  Z 

■<- 6  V,,  F  •  4*  ^  +  6  F^&’V,  4  +  £l7»4'V*& 

+  £  *eJ/tov,i.V* J,.  +txv,$”7h& 

*  «'/<•  V*  ‘H*  vh  5  +  *  Vl  V*  f,"  ■  <7*  £ 

+  -  >(%/„)!  ^  }&$)  l 

vj 

+  *3(^, £.  4" )  ?  +  «’/*.  fa  f„ •  v,  £ )  £ 
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+  ^  fa*  ju  ' ^ k  )  <k  *  £3  P fu  fa  ^  4 ) 
■*■  eV£*  fo  •**+«. )  +  (4a  )*  ■*•  &z^zJi 

Z  '  2- 

*•’££  * ♦  «**..££ 


r*  *? 


*  &il\  i«.s  <j>,  4 


o 


where 


F  *  /*•  "  *  jl;  c*j  /(i 


(?■!■“ i) 


The  expressions  ,  (V^) ,  and  gz  must  now  be  integrated  from 
the  bottom  to  the  instantaneous  free  surface  1^  .  The 
integrals  of  the  known  cross-space  dependencies  are 
expressed  in  terms  of  the  symbolic  values  G  and  A.  ,  which 
are  defined  in  Appendix  D  and  evaluated  where  neccessary. 
The  integration  results  in 

A 


i.  */i*A 
4. 


-  <5^ 

(bf-  A,L)  j.  t/7;  -4-  62<J 

1+W}*  (c./jp) 


'V  "v-* 


»s j* j'* j  'ji\'  V  .'l ,  ■■  . 


-hk 


■*  6i„+y )  4t 

“  (X  *£  )  ^  K>  *  ^  1  ^3.  ) 


£././* 


I  &£f  d* 
\  x 


-  6l^w  £  +e'l F'diCvJ,)1 
x  ■*  4  x 

+  ^(K+yXyuh' )*  +  £v^z  (^4)V  +  /V„4^  h"de 
z  x 

*  6  ^.4  •  [v*  4,  * 6  {&*  &j  ]v*4om  vh  & 

-V  J  ' 

A  0  ^ 

+  £  /  /'®  ^A  4  4  *  fel^lo  ^  <k  ’^i»  4 

-A. 


♦  &HvHh"'Vhh/)(f’.+rl)  *  t'A  vh£t"'rh  4 

+  *lJ)v-vh  4'  4  *  ^  4/-  ^  £ 

+  £f{^  £°*}  &o*j  «j  ^ 1  ‘ 

•  v*  £  ]  £  *■ 
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+  £*{&  -x£o/j 4j]rh$t-rM£  *  i'AvHb'-Vki. 
+  I [k"  )l  di  +  6l[ (Fil'd*  jjj,  +  6*  Jl2* 

^*4  X  —A  "  ^ 


/V  /-v 

+  *J  h 4>"z  d*  <k  *  t'jhiKz  <t*  k 

~l\  ~"^e 

+  ‘■1{&  -  >  .ff  “5  A*  1  $  ?» 


(ci.  a) 


Here  ,  the  remaining  integrals  are  given  by 

,n 


r. 


~h> 


d*  *  G>to  -  Za  £  £>». 

,\*i  ’>*  J 


1  3 


/*t 

J  d* 


3 

£  - 

j*» 


IvhF/i 


A*’6*; 

(c .  i.  a) 

(C.  (.<•<') 

7 

sj  ''i^°/jv^ 

CC.  /.  AT) 

expanded  about  the 

0(1)  mean 

%. : 


3; 

5 


i 


>3 


cd 


.’-  - i-. 


water  level  according  to  the  expansion 


G  a  G  ■+  £  i G  +  £1^(1G  +(iji+  bi,)  J  + 

+  ^  (*)i+bx)  &  'J  +  .. .  (c.i.k. 


as  shown  in  (D.l-3) ,  and  similarly  for  /A  .  The  integrals 
containing  ^  are  evaluated  here. 


*  j^(Vi V  ^ + £  fa  v„  )(fa  •  £ ))  * 

+  £\^KOv<^K  ft)  (c.i.n) 


where 


# 


«■  6 1  )) {- Vh-fyvh^' )]  &l  f 


hither  order  term 


(C.Ut) 


where 
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A. 


dl  *  i&rtn  it*  L  »f  0  ,  £  ) 


<vXv*)fv(^Ol 


holier  ctrlv 


Cc.1.  rt) 


vw 


l  i»\h  iz  *  w  -  )]lf,  +  AtyUr  Old*  ten 

(C.I.2o) 


-i. 


SH'*  *V„<fc"j  <?, 


whose 


5 


is 


•  % 


S 


-j 

w 


.  ^ 


•a 

s* 


tiBBeftgaflasass^^ 
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ft 


H 

= )  h> 

-ft 

We  also  note 

that  the  terms  ^Tej. 

and  Gowvi  may 

be  combined 

to  yield 

g;> 

+  C;  ■ 

*  f*ful  Vi.^ 

fc./.  21 ) 

The  Lagrangian  is  constructed  by  adding  the  terms  ‘ 
through  1^  and  expanding  the  integrals  G  and  Jk  .  The 
final  form  of  L  is  written  in  the  expanded  form 


where 


L  -  L.  *  &  L, 


+  6.'*  Ly  +  ... 


+  -  Y»  +&£?  ]  *  Q  (' n  •  (h  ^  4'))z 

- «?-  v»  ‘f,>*  fa  ■  (i  v>  )i 


+ 


“  %  %.  +  3  46^ +*>>.)  +  *  A  'Pjft 

+  fc  {  4't  *  4'-  r„  4'  -  4]  *  J"1 

+f&w  ■  %  4 )  4,  *  4t  4  ? 

4i  +  (ix+^{2LJth2k&ll 

x 

+(WH  -Oivki')jg  -  gs\  .  (wu  o  7,  l 

-  ty  *  4)  (vh&-vh)[vk-(£vh  4 )  J 

-(%&'•**  )fa-0>q;  4.')}  y,1 

■*■[<»•. _  ^ 4^»j  -JTI  °<i  V  ] fa4 ) 

<-  fe  -  2x£  -<j  £j  -  £  £  «j  <*fc  Gj*  J  J_ 

j**  hi  j*»  g  J  ^7 

-^Vh<f>,'  |Vh  &,j )  +  ^  /.j  J  ^ 


=  *  7,  { <K'4  ♦  v,  4.'-*  4J  - 

+  +  64-^<p/? 

+  {&„'->£ ^j^/]  -v*  <j>L  +  &Z -vkti  ft?, 

j*« 

-fX. •*»£)$  *  (X. -^4 )  ?«. 

+{c.:  -  6.J  ]  1,  (SiM  +(£, -V-J,  ){,£ 

>{£-zl^h£ 

+{£ -+k^J>-]vhh-V.h 

>W -+£**}'}  4  4 


+  or!1) 
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a  <j  +  ^4  -  YlJ  + 

■*“  ^  fe  bt  )  +  f2?|  +  ft  &f  )  * 

2. 

+  (Y7*.*fe>1.)/6  *ft^  ) §t^  +  7* 

+  ((*li+K)c,l>  ■nf.^rM)  +*itJfvk$,-rk& 
+{{yl  +  bx)Gim  +  ft&l')  ^  +  ^  W 

**'<k'*&  ^zw^w/r 


These  are  also  given  as  (3.2. 1-5)  in  the  main  body 
of  the  text.  Here,  the  operator  D/Dt  is  a  total  derivative 
following  the  0(1)  mean  flow  -  vk<  <*/♦  4 ) . 


C.2  The  Linear  Equation  for  $>  to  0(  £  ) 

Varying  L x  with  respect  to  <£>  and  performing  the 
partial  integration  leads  to  the  expression 

-&L  -fo-X)?,  *  Gw  I,  *  *  VfC.  ]  ?, 

*{&l  -  ut  £i  - kt  -i  6jk }  ?, 

Expanding  slightly  and  performing  some  cancellations  leads 
to  the  form 

~ Pfl.  f  fit  ?,  - %•  i  +  £6^,. 4i 

+  f2^~j  +  "<i  <*»  1  ^  ^  + 

*  <£j  ~  ZZ* j  ^  <ij>  J  <k  -O  (c.Z 2' 
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-k  ‘w.  n. 


■  y 

rt 

■  \ 

'  ':£ 


a 


•AS 

’W 


fl 


The  terms  involving  ,rv-j  are  0(S  ),  and,  in  those  terms, 

t<7> 

^ ^  may  be  approximated  by  -k  Cp>  .  Defining  the  quantity 

*  ^6)*  *  £j* 

we  rewrite  (C.2.2)  as 

-J2g  -fa-n'i'i,  *  ?,  -  v*  •£&.'.  ft?,?  * 

+  {^w. 4>,  +  2^ v; ( £ ■f’,  ♦ 

^  g 

v;|,  -  J  -  ££  i-o 


ng  to 


(C  2.3) 

r?  z'?’ 

The  term  /j,  ^  may  be  expanded  according  to 

Vfc’J  .  -fel£  +  2;frVH%  *o(sl)  (c. 2.v) 

where is  defined  in  section  (2.3)  and  is  related  to  the 
spatial  variation  of  the  amplitude  of  <£>  .  Noting  from 


Appendix  D  that 

J-s1  .  /•*'  - 


k  &oj  +  &o\  S  &  *  j 


i a V 


fc.2.>0 


it  is  apparent  that  all  0($  )  contributions  other  than  those 
provided  by  the  first  four  terms  of  (C.2.3)  dissappear. 


Also,  the  term 


is  0(S  );  Vh  can  be  approximated  by  .  The  term 


t 

VoVt 


is  identical  to  the  higher  order  term  found  by  Smith  and 
Spr inks (1975) ,  after  absorbing  a  boundary  term  through 
Leibnitz  integration.  Using  (C.2.4)  and  (C.2.5),  (C.2.3) 
can  be  written  in  the  form 


“fe  '  ~  h  <j>( 

+  ^  l  jT,  ^0  }  fti 


=  O 

(c.2.0 


Equation  (C.2.6)  is  the  correct  equation  for  S  to 

Tl 

0(£  );  however,  the  operator  is  only  symbolic  and  is  not 

/V 

practically  applicable.  Letting  <j?(  be  given  by 


AP'jw***: 


V  \* 


3 


we  can  write  in  the  form 

Vh£  •  -^-(vHA)e 


(¥)« 


(c.2.7) 


where  Vj,A./A.  is  one  component  of  o^.  In  order  to  make 
further  explicit  progress,  it  is  neccessary  to  expand 
into  its  components,  and  then  to  take  the  derivatives  of  the 
products  c-ij  .  We  define  the  quantities 

■*  Ct^)  j  L  »  fe-V|.  C  ■ 

k  2kc 


/Ss  =  ;  &  *  is  vi.  (&*•*')  ■ 

k'A  ftV  ’ 


*  k-yfa-’ll) 

1 

where  -  £>$■  are  0(  S  ),  and 

&t=  j 


r  fe' 


£7=  1st 
feV 


(c.  2. 0 


6.  z.O 


</  «*  */  <  *•*  •  *  k  •  «  *  •  *  .*  <  ’  •  *  «*  .  ‘  ,*  .*  ,»  •  .* 


:a 

s  a 

J 

: 


where  j^t  are  °(  >  )•  Also,  we  let 


(C.LIo) 

The 

coefficient  c2. 

L  may  then  be  expanded 

as 

^ i  * 

4* 

~  ~  ^3 

(C.  2.U  ) 

The 

terms  may  be  written  as 

J?  •*»<*,  *  - 

2/?  ex',  0^3 

+  *>, 

• 

1-%*X  -  - 

2fe2c^3  + 

(^3  ~^y 

i-fg- 

r/4)  - 

— 

y>V»^  + 

Hkl<^3jS7  + 

Ck*o<: 

J  +  k  f 7 

£’^3  *  * 

~ok  Vj  +• 

) 

(c.  2,li ) 

and 

the  terms  k»  vl 

G*j  33 

55 

(c,  2.  n) 

where 

=  —  O’  kh  kit  (l  -  fell  7/**/  kb  )  ? 

^  fimh  2MM  C  J 

2a,  =  fCZ-ty'K  )  7 lahL  kb  +  2kh  [%  -£•)[/  “  2  (/-  &, 


where 


a 

-cr 

ej  S/*A  Zkh 


fe  ’  ^03  "  +  °^?  "3*4  C^‘  2'  l' 

where 

3*.  .  *(z*-i)-kJ,kh  7 

^  ^  3  S/U2^  J 
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*.«.  =  - 


0  *  J&uLkk  4- 
^  1  SWJflM 

■f  f  (hit)  _ -  f  /-  kh  IahI  kh')~\ 

3  jvi>i  2-^  J 


Substituting  (C. 2. 7-15)  into  (C.2.6)  leads  to  the 


equation 


+  ^“^A?  +/^y  A  “  °^i 

+  ^2  <*3*  ^  -ft  )  °4  +  (~^ot  ^  4/ 

+  +  ^ZZ  (fr 'ft  ~  fo  H )  =  ° 


(c.  2,  IL  ) 


where 


< n.s.y  * 


'S' 


/V  / 

-V‘h-(c>fVh4>i  )  +  4  (c.2.n') 


'W 


V 
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and 

y,  -  9„  + J 

£  '  -  3.,  -  ?.r  '  3„) 

y3  *(*2,  -2t-2„  -22*- Vr) 
X,  *  3,,  ♦  >„ 
yr  •  V21  -  2  3* 

4  ’  3U  ~  3«.  *  3" 

and  where  we  have  used  the  results 

5|l  +  *  ^21  +  *  ^23  ^  ^3  *  ° 

Note  that 

(  n.s.y  =  o 


Cc.j.ff) 


is  equivalent  to 
dropping  a  small 


Booij’s  mild  slope  wave  equation  after 
%z 

0(  i  )  term  which  appears  in  his  equation 


ta 


at 


due  to  an  ertoc  in  the  derivation.  After  substituting  the 
free  surface  boundary  condition  (3.3.7)  into  (C.2.16),  we 
obtain  the  reduced  wave  equation  for  ,  correct  to  0(  5  ), 

given  by  (3.3.8),  where  an  arbitrary  dissipation  term  is 
included.  The  coefficients  2f|  -  ^  are  given  explicitly  in 
Appendix  D. 


APPENDIX  D.  INTEGRALS  APPEARING  IN  L 


In  the  derivation  of  the  Lagrangian  L,  various 
integrals,  denoted  by  G  and  Jl.  ,  arise  in  the  integration 
over  the  depth.  The  notations  used  are  as  follows. 

Integrals  of  fj  quantities: 


etc . 


Integrals  of  fJfl  : 


Integrals  of  products  of  f  and  f ^  quantities: 
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**  ’  4*  -  /  {» %  f.* 


*  ^tc • 


Bach  of  the  integrals  G  and/fi  may  be  expanded  about  the  mean 
water  level  z*ba  by  the  following  procedure 

?  fc. 

X- 1  Ad*  *  Ad*  +  Ate 

“A*  **4b 

6  6[t  +  (;x(jt+ki) 

*  I  Ate'  f  /  A  rf*'  (p,i) 


where  A  is  an  arbitrary  function,  and  the  wave  reference 


level  z  has  been  shifted  to  the  slowly  varying  0(1)  mean 


water  level.  The  coefficients  f ^  ,  fZe  are  then  defined 


with  respect  to  the  local  mean  depth  h=hft+btf  rather  than  the 
still  water  depth  h6  .  The  integral  from  z»b6  (z'*0)  up  to 


the  local  water  surface  is  further  expanded  to  give 


X  9  Ad?'  4*  A i  +  fa,  *6Yy,+s>a)]  A^i  t 


*£Mi±  *Y/>.+OS  Aaa  j  +  " 


a  I 


*  * 


4 


-I'*  U?.  +  « 

ftu) 


The  contributions  to  I  are  then  reordered  in  powers 
give 


r  -  i'  *651,1'  +  4l{^1+fe,)x"*  i,'x 

*  6!{2?,^*t,)I**  fl1!®]  *  ... 


/// 


of  £  to 

)* 


4>j) 


The  integrals  appearing  explicitly  in  the  equations  follow. 
For  the  mild  slope  approximation  (terms  of  0(  S  )  in  «jJ, 
neglected),  integrals  of  the  form  G0 ,  G0#  ,  and  A  appear 
in  and  L^.  These  are  given  by 

*  1 

&l>mCC\h  <£l  =  el  Cl- *■)/<& 

where  n  is  the  ratio  C^,  /C  given  by 

'H.  =  4?  (  l  +  }  5 


8 


kt 


~  2  k  CosL  kk 

•S' /*/  *kh 


®  1 _  i"  C  ox  It  Zkh  +  cos  l>  kl  1 

Z  cosh  kh  S/*h  ^kh  C  1 


a." 

A  " 


Coil  kit  S»1H kh 


COsL  Zklt  -  cUkl 


1  > 
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The  remaining  required  integrals  are  the  quantities  G 
i,j«l,2,3,  and  the  combinations  3^'  given  by 


which  are  given  by: 
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The  combinations 
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are  equal  to  zero.  Referring  to  the  definition  of  F 
(3.3.8),  given  by  (3.3.10),  and  to  (C. 2. 13-16),  t 
remaining  coefficients  are  given  by 
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The  coefficient  U  ,  expanded  in  (3.3.13),  has  the 
terms  ^ ,  which  are  given  by 
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APPENDIX  E.  CONSERVATION  EQUATION  FOR  WAVE  ACTION 


In  order  to  derive  several  of  the  results  in  Chapter 
4,  an  equation  governing  the  flux  of  an  energy  related 
quantity  is  required  in  addition  to  the  Euler  equations 
found  in  Chapter  3.  It  has  been  shown  that,  for  propagation 
in  a  moving  medium,  the  appropriate  conserved  quantity  is 
the  wave  action,  given  by 

WAve  Acbo/u  »  B / <r  (El) 

(Garrett (1967) ,  Bretherton  and  Garrett (1968 )  ),  where  E  is 
given  by  E»g  JA[  /2 .  Note  that  the  expressions  for  conserved 
quantities  are  formulated  here  with  respect  to  unit  masses, 
which  is  possible  since  the  density  has  been  assumed  to  be 
constant  throughout  the  fluid  domain.  The  lowest  order 
conservation  equation  for  E/d-  may  be  obtained  directly  from 
the  linear  wave  equation  (3.4.1),  using  the  definition  of  W 
given  by  (4.2.1): 

#  +fa- it)  &£  -v*-(cc jV*  ?,)  +  ?'  " 

=  a. 
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where  w  in  the  last  term  is  a  coefficient  related  to  the 
rate  of  energy  dissipation  (  Booij(1981),  Dalrymple,  Kirby 
and  Hwang  (1983)).  Following  the  approach  of  Jonsson (1981) , 
we  assume  that  <p,  may  be  represented  by  the  ansatz 

X  <j>,  *  &  (£■’$) 

where  R  and  9  are  real  quantities.  Comparing  (E.3)  to 
(A. 4),  it  is  apparent  that 

R-^-jAl  (e«) 

while  $  absorbs  the  higher  order  phase  effects  neglected  in 
the  definition  of  *Y'  .  Since  we  are  looking  for  the  lowest 
order  results,  we  may  approximate  9  and  its  derivatives 
according  to 

9  -  OC (F.Sk) 

=*  i lt  +  o  c 6l)  ~  (e.^0 

vhe  *  vMf  +  o(6l)  =  j?  (£.?*) 

Substituting  (E.4)  and  (E.5)  into  (E.2)  leads  to  the  complex 
equation 

+  (yk-li)  S£  -  vr  (ccyvkR)  (<rtR  +  l<rRt) 

--if  Vk{!£r)R  *-2(r!i)-VkR  ■k%'(J'Cc>)R  -  2£cc}  •  CJ  R j 
-yi  C~wR  m  O  (£'£) 
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Since  R  and  Q  are  real  by  definition,  the  real  and  imaginary 
parts  of  (E.6)  must  individually  equal  zero.  Gathering  the 
imaginary  terms  and  multiplying  by  R,  we  obtain  the  equation 

-o-u-vi,/?'  -vk-(o-U)  Rl  -cc^t-v,Rl 


(kccy~)  Rl  -  cwR1  -o  (f-7) 

Then,  noting  that  we  can  write 

cc^Jz  -  kc(c^k/k)  =  crQ^  (e.9) 

we  rearrange  (E.7)  to  the  form 

C<rR\  +  V(,  •  (<rR'(a+&>))  *  -<rwRl  (e.i) 

* 

Now,  R  may  be  replaced  by 

R1*  cfiA\X/<rl  «  -&(■-)  (E.io) 

which  yields  the  relation 

(e/r\  i-  =  -w(eA)  (< r.ii ) 

It  is  apparent  that  w>0  for  a  dissipative  motion. 


Neglecting  the  right  hand  side  of  (E.ll)  gives  the  standard 
equation  for  conservation  of  wave  action  in  a  conservative, 
inhomogeneous,  moving  system,  given  originally  by 


Garrett (1967) .  The  result  with  dissipation  has  been  shown 
to  be  somewhat  incomplete  by  Christoff ersen  and 
Jonsson (1980)  ,  who  show  that  the  right  hand  side  should 
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contain  a  term  involving  the  bottom  shear  stress  resulting 
from  the  0{1)  mean  current.  In  this  study,  we  regard  the 
ambient  current  as  an  imposed  flow  and  neglect  frictional 
effects  related  to  the  current,  thus  yielding  the  version  of 
the  conservation  equation  given  by  (E.ll) 

Jonsson (1981)  shows  how  the  equations  governing  the 

geometric  optics  approximation  of  (3.4.1)  may  be  obtained  by 

using  the  real  part  of  (E.6).  In  this  case,  it  is  necessary 

J. 

to  retain  the  expressions^  ,  and  f.  in  order  to  obtain  the 


*  desired  results. 


APPENDIX  F.  BOUNDARY  INTEGRAL  METHOD  FOR  THE 
ONE-DIMENSIONAL  PROBLEM 

The  one-dimensional  problem  (topographic  variation 
in  the  x-direction  only)  treated  in  Chapter  5  admits  an 
"exact"  solution  by  means  of  the  boundary  integral  method 
(BIM) .  The  method  is  exact  in  terms  of  its  basic 
formulation,  whereas  the  theory  of  the  main  body  of  the  text 
relies  on  series  expansion  representations  for  the 
propagating  components  of  the  wave  field  alone.  However, 
solution  by  BIM  still  requires  discretization  of  the  domain 
and  numerical  evaluation  of  the  results. 

The  method  used  in  this  study  was  described 
originally  by  Lee(1971)  and  was  applied  by  Raichlen  and 
Lee(1973)  to  the  study  of  linear  waves  radiated  in  one 
direction  from  an  infinite  inclined  plate  wave  generator. 
Here,  we  extent  the  method  to  waves  travelling  at  an  oblique 
angle  to  the  direction  x.  Let  the  total  velocity  potential 
^  be  represented  by 


jlL* K.^.  -  Uit) 


(F.l) 


where  m*0  for  normal  incidence  and  0<m<k^  for  oblique 
incidence,  where  is  the  wavenumber  for  a  constant  depth 
region  A  on  the  upwave  (x>  -  oo  )  side  of  the  topographic 
obstacle.  The  incident  wave  ^  is  then  given  by 


<Pz  ’  £x  £ 


(Fi) 


with  unit  amplitude,  where 
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Substituting  <j>  into  Laplace's  equation  (2.1.1)  leads  to 


modified  Helmholtz  equation 
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with  the  associated  free  space  Green's  functions 
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where  K_  is  the  modified  Bessel  function  of  the  second  kind 

O 

and  zeroth  order  and  r  is  the  radial  distance  from  the 
singularity  to  the  field  point  in  question 


f  1  t  7  K- 

r  *  I  (x-Xt,)  *  (z-- O  ] 


A  fluid  domain  Ji. (x,z)  ,  shown  in  Figure  F.l,  is  constructed 
by  separating  the  constant  depth  regions  from  the  region 
containing  the  variable  topography,  but  including  enough  of 
the  constant  depth  region  to  either  side  of  the  variation  to 
allow  for  the  decay  of  nonpropagating  modes.  The  solutions 
in  the  exterior  regions  A  and  C  may  then  be  given  by 


it ^  “  /. ^  R  (cocal*  Coih  fynfa+d)  <£ 


Cm* 


It  •  -^T  (  CO  COi  L  klc)  call  kcChc+z)  e  (F,7lo 


where  R  and  T  represent  reflection  and  transmission 
coefficients,  and  are  in  general  complex  to  allow  for  phase 
shifts  in  x  and  y.  Conservation  of  energy  in  the  diffracted 


Figure  F.l  Boundary  integral  method  domain 

in  relation  to  extent  of  variable 
topography 


wave  field  requires  that 


1*1  *  hi  - 


Also,  the  boundary  of  J2  is  denoted  by  oSl  •  Following  the 
standard  procedure  for  constructing  solutions  by  means  of 
Green's  functions,  we  take  the  inner  product  to  be  given  by 


v*.  il  *  il 


(F.O 


Then,  using  Green's  Second  Identity,  a  formula  for  the  value 
of  3l  (xd,zft  )  for  an  interior  point  of  the  domain  fL  is 


obtained: 


I(\ o -  f  fl© &GHS0  -  (Fte) 

an. 

where  n  denotes  the  outward  normal  to^xi.  Values  of  in 
the  interior  of  ai_  are  thus  completely  determined  by  31  on 
>11,  and  the  problem  is  reduced  to  finding  31  (  iJ"_)  .  Values 


)  to 


of  5  on  J 0.  are  determined  by  moving  the  point  (xd,i6 
the  boundary  and  determining  the  principle  value  of  the 
integral  in  (F.10)  (see  Figure  F.2) .  Including  the  point 
(x. »z„)  within  the  contour  and  retaining  the  residue  leads 

o  o 

to  the  formula 

I Cx.)  *  if  - &cag(x)  ] jLsct) 

ill 

for  points  (x&  ,zfr)  on  SCL  .  Substituting  for  G  from  (F.5) 
leads  to  the  formulas 

2  60  =  {  F.(*r)  * (X)  -  2CZ)±(K(~))1  dSlx) 

id  (F.lU) 

for  oblique  incidence/  and 


U'i.)*-±.f{Ar  Sex)  -ZC&' ll.fartfjsb) 


(f.iil) 


for  normal  incidence. 


The  details  for  the  numerical  approximation  of 
(F.lla)  follow;  the  corresponding  details  for  (F.llb)  may  be 
derived  in  analogous  fashion  or  obtained  from  Raichlen  and 
Lee(1978)  .  First,  the  boundary  is  discretised  into  N 
boundary  segments  which  are  presumed  to  be  straight  line 
segments  (thus  leading  to  discretization  errors  in  the  case 


of  locally  curved  boundaries).  Taking  x^  a{x0'z4J  t0  be  the 
boundary  point  for  which  _2  is  being  determined  and  x j  to 
vary  over  all  boundary  points,  (F.lla)  is  approximated  by 

j  (&) =tt<L  .{  k  (*o,  )j£cti>-i.aii)kLz  (■*** ))? 

JSI  (£.12. 

where  c.s*.  is  the  length  of  the  jtn  boundary  segment  and  .£ 
is  assumed  to  be  uniform  over  each  boundary  segment. 
Special  treatment  is  needed  for  i  =  j,  where  r^  *0.  Writing 
(F.12)  out  in  matrix  form,  we  get: 
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and  the  evaluation  of  (F.14)  for  i^j  is  straightforward. 
For  i=j,  G  .  •  in  a  small  neighborhood  of  x  •  can  be 

"A 

approximated  by  (Abromowitz  and  Stegun  (1972)  ,  p.375) 

<iSx  CF.< 

Integrating  over  the  boundary  segment  leads  to 


<■  v  1.  '  J 

Following  Lee(1971)  and  Raichlen  and  Lee(1978) ,  G„  . .  is 

*  *+ 

given  by 

(P.  n) 

due  to  the  neglect  of  the  local  curvature  of  the  boundary. 
\*»  A 

Finally,  •*-!;  can  be  specified  in  terms  of  boundary  and 
radiation  conditions: 
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leading  to  the  final  matrix  equation 
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where  I  is  the  diagonal  identity  matrix 


Kirby  and 

Dalrymple ( 1983a)  have  compared  the  results  of  the  present 
approximation  to  results  obtained  using  integral  equations 
based  on  eigenfunction  expansions  of  J[  in  regions  of 
constant  depth  separated  by  vertical  boundaries,  as  shown  in 


V 


2 
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Figure  F.3.  Two  such  comparisons  are  shown  in  Figures  F.4 
and  F.5.  In  Figure  F.4/  results  are  for  the  case  of  an  .  ; 

asymmetric  trench,  with  the  trench  being  twice  as  deep  as 
the  incident  wave  region,  and  the  transmitted  wave  region 
being  half  as  deep  as  the  incident  wave  region.  The  ratio  r* 

of  trench  width  to  incident  region  water  depth  is  5:1,  and 
only  normally  incident  waves  (m»0)  are  studied.  Generally 
good  agreement  is  found  between  the  BIM  and  eigenfunction 
solutions.  The  results  of  Lassiter (1972) ,  included  for 

comparison,  are  seen  to  deviate  significantly  from  the 
solutions  of  Kirby  and  Dalrymple (1983a) .  In  Figure  F.5,  BIM 
results  are  shown  for  a  symmetric  trench  and  an  oblique 
angle  of  incidence  of  45^.  Here,  h^/hj =3  and  L/h,*10,  using 
the  notation  of  Figure  F.3.  For  values  of  kh, <0.792  for  the 
incident  waves,  the  value  of  the  x-component  of  the 
wavenumber  vector  in  the  trench  is  imaginary;  this  result  is 
seen  to  present  no  difficulties  to  the  analysis  technique. 


The  method  used  in  the  present  study  represents  a 
relatively  crude  approach  to  the  boundary  discretization. 
For  a  more  satisfactory  approach,  where  Si  is  assumed  to 
vary  linearly  over  each  segment,  see  Liu  and  Abbaspour 
(1982)  . 
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Figure  F-4. 


Reflection  coefficient:  asymmetric 
trench 

hj/h^  »  2,  L/h^  -  5,  *  0.5 

•  x  -  0. 

K-  -  |R| ,  K  -  w2/g 
*  Kirby  i  Dalrymple  (1983a) ; 

- , -  Lassiter  (1972); 

•  BIN  solution 

(from  Kirbv  and  Dalrymple,  1983a) 


Figure  F.5  Transmission  coefficient:  symmetric 
trench 

hj/hi  *  3,  L/h^  =  10 


-  8^  »  0°;  -  8^  *  45° ,  Kirby  and 

Dalrymple  (1983a) 

<  BIM  solution 

(after  Kirby  and  Dalrymple,  1983a) 


APPENDIX  G.  RELATIONS  FOR  THE  DISSIPATION  COEFFICIENT 
w  BASED  ON  VARIOUS  DISSIPATION  MECHANISMS. 


In  this  appendix,  the  cole  of  the  dissipation 
coefficient  w  in  the  governing  equation  for  wave  motion  is 
investigated,  and  its  value  is  related  to  physical 
quantities  based  on  various  mechanisms  for  wave  energy 
dissipation. 


The  role  of  w  in  determining  the  wave  amplitude  can 
be  investigated  by  using  the  current-f ree  form  of  the 
elliptic  model  (4.2.4)  ,  wnich  reduces  to 

Vh'  (CC%  4>,  *  O  C&-I) 


Considering  the  simple  case  of  waves  propagating  in  the 
x-direction  over  a  flat  bottom  reduces  (G.l)  to 


:: 


where 
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kc  s  k  (  /  +  )  6S.3) 

The  solution  of  (G.2)  for  waves  travelling  in  the  positive  x 
direction  is  given  by 

/%»  vc  )C 

<M) 

where  A  is  real  and  arbitrary  or  can  be  related  to  an 
initial  condition,  if  kc  is  expressed  in  real  and  imaginary 
parts,  kc  a  kr  +^k^ ,  the  exponentially  decaying  solution  for 
is  given  by 

/v  —  fejC  X  J.  &r-  X 

(j>(  Cx)  »  A  e  c  &.f) 


Solving  (G.6b)  for  w  yields 
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Depending  on  the  nature  of  the  localized  energy 
dissipation,  w  takes  on  different  values.  For  a  number  of 
dissipative  mechanisms,  the  wave  height  decay  is 
exponential,  in  agreement  with  (G.5) ,  and  k^  is  known. 
Several  examples  are: 


1.  Porous  bottom  (Reid  and  Kajiura,  1957). 


2kk 

•+  still*  2kh 


fax) 


2.  Viscous  mud  bottom  of  thickness  d,  viscosity 
,  and  density  ^  (Dalrymple  and  Liu,  1978) 


M'Li'imn 


I'-rMt'-  f  (&. 

where  k;  is  for  shallow  water  (for  simplicity) 
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3.  Laminar  bottom  boundary  layer  (Ippen,  1966) . 

k:  *  kZ  (&A0) 

/fl  ZkA 

4.  Densely  packed  surface  film  (Phillips,  1977) . 

,  .  Vt 

x  ( v/lu)  - 

n  2  7^4  kA 

For  case  1,  a  more  complete  expression  has  been  obtained, 
including  the  effect  of  boundary  layers,  by  Liu(1973). 

For  other  types  of  damping,  the  wave  amplitude  is 
not  exponential  but  varies  as 

A  60  *  Ac  _  fc.n) 

l  +  <*x 

where  is  the  initial  wave  amplitude  and  si  is  a  damping 
coefficient.  For  a  rough  bottom  with  a  turbulent  boundary 
layer,  Madsen(1976)  showed  that 

ci  9  If*  _ k  At> 

3tt  s' mb  kb  (2Jeb  i*”  s>*h  2^ A') 

where  f^  is  a  wave  friction  factor.  Dalrymple,  Kirby,  and 


Hwang (1983)  derived  a  dissipation  formula  for  waves 
propagating  through  an  array  of  vertical  cylinders  occupying 
an  arbitrary  fraction  of  the  total  water  depth,  based  on 
drag-dominated  dissipation  for  each  cylinder  times  the 
density  of  cylinders  per  unit  area.  The  resulting  value  of 
(X  is  given  by 

where  s  is  the  position  of  the  top  of  the  cylinder  relative 
to  the  bottom,  b  is  the  spacing  between  cylinders,  D  is  the 
cylinder  diameter,  and  Cp  is  a  representative  drag 
coefficient  for  each  cylinder  acting  alone. 

In  order  to  relate  x  to  w  it  is  convenient  to  first 
relate  <~X  to  k^  by  comparing  the  following  expansions 

+  *  ]-o^X+  (olx)1  - 

C  *=  I  -  kx  +  fax)  (kxx)\ ...  C&.tfb) 

1  "  3  " 

Therefore,  for  small  values  of  k^x  or  x,  we  may  take 
can  thus  use  (G.7)  in  either  case  with  c. 


±  IsMs  )  ^ 

ZsijJi I  kh  ( 2  kh  +  siih  J 


replacing  ,  as  long  as  oi  Ax  is  relatively  small.  For 
subsequent  numerical  modeling  it  is  sufficient  to 
approximate  (G.7)  by 
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